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Abstract 


The high-charge-and-energy (HZE) transport computer program 
HZETRN is developed to address the problems of free-space radiation 
transport and shielding. The HZETRN program is intended specifically 
for the design engineer who is interested in obtaining fast and accurate 
dosimetric information for the design and construction of space mod- 
ules and devices. The program is based on a one-dimensional space- 
marching formulation of the Boltzmann transport equation with a 
straight-ahead approximation. The effect of the long-range Coulomb 
force and electron interaction is treated as a continuous slowing-down 
process. Atomic ( electronic ) stopping power coefficients with energies 
above a few A MeV are calculated by using Bethe’s theory including 
B ra gg s rule, Ziegler 's shell corrections, and effective charge. Nuclear 
absorption cross sections are obtained from fits to quantum calcula- 
tions and total cross sections are obtained with a Ramsauer formalism. 

Nuclear fragmentation cross sections are calculated with a semiempir- 
ical abrasion-ablation fragmentation model. The relation of the final 
computer code to the Boltzmann equation is discussed in the context of 
simplifying assumptions. A detailed description of the flow of the com- 
puter code, input requirements, sample output, and compatibility 
requirements for non-VAX platforms are provided. 

1. Introduction 

During the last 40 years, propagation of galactic ions through extended matter and determination of 
the origin of these ions have been the subject of many studies. Peters (ref. 1 ) used the one-dimensional 
equilibrium solution, without including ionization energy loss and radioactive decay, to show that the 
light ions have their origin in the breakup of heavy particles. Davis (ref. 2) showed that one-dimensional 
propagation is simplistic and that leakage at the galactic boundary must be taken into account. Ginzburg 
and Syrovatskii (ref. 3) argued that the leakage can be approximated as a superposition of nonequilib- 
rium one-dimensional solutions. The solution to the steady-state equations was given as a Volterra 
equation by Gloeckler and Jokipii (ref. 4), which was solved to the first order in the fragmentation cross 
sections by ignoring energy loss. This provided an approximation of the first-order solution that 
included ionization energy loss and was only valid at relativistic energies. Lezniak (ref. 5) gave an over- 
view of the cosmic ray propagation and derived a Volterra equation that included the ionization energy 
loss and evaluated only the unperturbed term. The previous discussion indicates, that for a long time, 
the main interest of cosmic ray physicists was to achieve first-order solutions in the fragmentation cross 
sections where path lengths in the interstellar space are on the order of 3 to 4 g/cm 2 . Clearly, higher 
order terms cannot be ignored in accelerator or space shielding transport problems. (See refs. 6-9.) 
Besides this simplification, previous cosmic ray models have neglected the complicated three- 
dimensional nature of the fragmentation process. 

Several approaches to the solution of high-energy heavy ion propagation that include ioniza- 
tion energy loss have been developed (refs. 6-19) during the last 20 years. All but one (ref. 6) have 
assumed the straight-ahead approximation and velocity-conserving fragmentation interactions. Only 
two (refs. 6 and 9) have incorporated energy-dependent cross sections. The approach by Curtis, 
Doherty, and Wilkinson (ref. 14) for a primary ion beam represented the first-generation secondary 
fragments as a quadrature over the collision density of the primary beam. Allkofer and Heinrich 
(ref. 15) used an energy multigroup method in which an energy-independent fragmentation transport 
approximation was applied within each energy group after which the energy group boundaries were 
moved according to continuous slowing-down theory. Chatteijee, Tobias, and Lyman (ref. 16) solved 



the energy-independent fragment transport equation with primary collision density as a source and 
neglected higher order fragmentation. The primary source term extended only to the primary ion range 
from the boundary and the energy-independent transport solution was modified to account for the finite 
range of the secondary fragment ions. 

Wilson (ref. 7) derived an expression for the ion transport problem to the first-order (i.e., first- 
collision) term and gave an analytical solution for the depth-dose relationship. The more common 
approximations used in solving the heavy ion transport problem were examined further by Wilson. (See 
ref. 6.) The effect of conservation of velocity on fragmentation and on the straight-ahead approximation 
was found to be negligible for cosmic ray applications. Solution methods for representation of the 
energy-dependent nuclear cross sections were derived. (See ref. 6.) Letaw, Tsao, and Silberberg 
(ref. 17) approximated the energy loss term and the ion spectra by simple forms for which energy deriv- 
atives were evaluated explicitly. The resulting ordinary differential equations in terms of position were 
solved analytically. This approximation results in the decoupling of motion in space and a change in 
energy. In Letaw’s formalism, the energy shifts were replaced by an effective attenuation factor. Wilson 
(ref. 8) added the next higher order (i.e., second-collision) term. This term was found to be very impor- 
tant in describing 20 Ne beams at 610 A MeV. The three-term expansion was modified to include the 
effect of energy variation of the nuclear cross sections. (See ref. 9.) The integral form of the transport 
equation was also used to derive a numerical marching procedure to solve the cosmic ray transport 
problem. (See refs. 6 and 12.) This method could accommodate the energy-dependent nuclear cross sec- 
tions within the numerical procedure. Comparison of the numerical procedure with an analytical solu- 
tion of a simplified problem (refs. 12 and 13) validated the solution technique to approximately 
1 -percent accuracy. Several solution techniques and analytical methods have also been developed for 
testing future numerical solutions of the transport equation. (See refs. 18 and 19.) More recently, an ana- 
lytical solution for the laboratory ion beam transport problem has been derived with a straight-ahead 
approximation, velocity conservation at the interaction site, and energy-dependent nuclear cross 
sections. (See ref. 10.) 

From an overview of these past developments, the applications are divided into two categories: a 
single-ion species with a single energy at the boundary and a broad host of elemental types with a broad 
continuous energy spectrum. Techniques, which will represent the spectrum over an array of energy 
values, require vast computer storage and computation speed to maintain sufficient energy resolution 
for the laboratory beam problem. In contrast, analytical methods (ref. 6), which are applied as a march- 
ing procedure (ref. 12) have similar energy resolution problems. This is a serious limitation because a 
final (i.e., production) high-charge-and-energy (HZE) computer code for cosmic ray shielding must be 
thoroughly validated by laboratory experiments; some hope exists of having a single code which can be 
validated in the laboratory. (See refs. 9, 10, 20, and 21.) More recently, a Green’s function has been 
derived which has promise for a code which can be tested in the laboratory and used on space radiation 
protection applications. (See ref. 22.) 

In this paper, the starting point is the derivation of the general Boltzmann equation. By using stan- 
dard assumptions to derive the straight-ahead equation in the continuous slowing-down approximation 
and the assumption that heavy projectile breakup conserves velocity, the Boltzmann equation is simpli- 
fied. A numerical procedure is derived with the coupling of heavy ions to the nucleon fields. Numerical 
stability and error propagation are discussed. The environmental model required as input to the HZE 
transport computer program HZETRN is briefly discussed. Atomic and nuclear models used to obtain 
the transport coefficients are discussed. Monte Carlo results are compared with the numerical proce- 
dures and database. Sample results for solar minimum and maximum periods are provided. Detailed 
descriptions of the flow of the computer code, input requirements, sample output, and compatibility 
requirements for non-VAX platforms are provided. This program, which is designated LAR-15225, is 
available through NASA’s software technology transfer center COSMIC (Computer Software Manage- 
ment and Information Center) at 382 E. Broad Street, University of Georgia, Athens, GA 30602. 
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2. Derivation of Boltzmann Equation 

Because the volume of any material is mostly electrons, the interaction of energetic ions passing 
through any material is primarily with these electrons. The cross section for the interactions of electrons 
is a at ~ 10“ 16 cm 2 . The long range of the nuclear Coulomb field also presents a sizable cross section of 
o c = I0- 19 cm 2 to the passing ion. Ion collisions are dominated by these two processes, but individual 
collisions have little effect on the passing ions. 

Although most collisions in the material are Coulomb collisions with orbital electrons and nuclei, 
the rare nuclear reactions are of importance because of the significant energy transferred in the reaction 
and the generation of new energetic particles. The transfer of kinetic energy into new secondary radia- 
tions occurs through several processes such as direct knockout of nuclear constituents, resonant excita- 
tion followed by particle emission, pair production, and possible coherent effects within the nucleus. 
Through these processes, a single particle incident on a shield may attenuate through energy transfer to 
electrons of the media or generate a multitude of secondaries which cause an increase in exposure. The 
process that dominates depends on energy, particle type, and material composition. 


The relevant transport equations arc derived on the basis of conservation principles by considering a 
region of space filled with matter described by appropriate atomic and nuclear cross sections. In 
figure 1 , a small portion of such a region enclosed by a sphere of radius 8 is shown. The number of par- 
ticles of type j leaving a surface element 5 2 dQ is given as <|>(.r + 8£2, £2, fjs 2 dQ, where 

£2, fj is the particle flux density, x is a vector to the center of the sphere, £2 is normal to the sur- 
face element, and E is the particle energy. The projection of the surface element through the sphere cen- 
ter to the opposite side of the sphere defines a flux tube through which pass a number of particles of 

given as x - 5£2, £2, E Js 2 di 2, which would equal the number leaving the opposite face if the 

tube defined by the projection were a vacuum. The two numbers of particles differ by the gains and 
losses created by atomic and nuclear collisions as follows: 

<|r(x + 8£2, £2, e]s 2 d£2 = £2, fjs 2 dQ 

+ 8 2 dQ^ £2', E, £')<!>*(* + K2, £2', E^jdQ dE 

k 

-5 2 dnfdl a. (E) + /£2, £2, e) {2A) 


where Oj (E) and £1 , £, E J are the media macroscopic cross sections. The cross section 

£2 . E, E J represents all those processes by which type k particles moving in direction £2' with 

energy E produce a type; particle in direction £2 with energy E. Note that there may be several reac- 
tions that could produce this result, and the appropriate cross sections of equation (2.1) are the inclusive 
ones. The second term on the right side of equation (2.1) is the source of secondary particles integrated 

over the total volume 28^ 8 2 d£l j y and the third term is the loss through nuclear reaction integrated over 

the same volume. The expansion of the terms of each side and retention of the terms to order 5 3 explic- 
itly result in 
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( 2 . 2 ) 


8 2 dQ [4l[x, £2, E ] + 8£2 • V^l, £2, eJ] = 5 2 dQ [fyfx, £2, eJ - 8£2 • V^[i, £2, eJ 

+ 26^ Jo. t ( £2, £2', E, £')^(i, £2 , E'jdQ' dE 

k 

-28 Gj(E)^x, £2, e)] +0(8 4 ) 

which can be divided by the cylindrical volume 28^ 5 2 d£2 J and written as 
£2 • Vfyfl, £2, eJ = ^Jct^£ 2, Q, E, E'Jtj)^, Q, E'Jdfl d£'-a.(E)(() y [l, £2, eJ + (?(8) (2.3) 

k 

for which the last term O (8) approaches zero in the limit as 8 -> 0 . Equation (2.3) is recognized as the 
time-independent form of the Boltzmann equation for a dual-species tenuous gas. Atomic collisions pre- 
serve the identity of the particle, and both terms on the right side of equation (2.3) contribute. The dif- 
ferential cross sections for the atomic processes have the approximate form 

2, £2 , E, E') = (£') S^Q • £2' - 1 Js^E + e„ - £') (2.4) 

n 

where n labels the electronic excitation levels and e n represents the corresponding excitation energies, 
which are small (1 to 100 eV) compared with the particle energy E. The atomic terms can then be 
written as 


£2', £, E'J<|> A (i £2', E')dQ dE -of (E)^[x, £2, £ J 
k . 

= X a y « {E + £ »> %{*’ ^' E + 0 ~ °f (£) */(*• £ J 

« 

- £a« <£> ♦/*. [oj: < £ > 5. e)] - of (*> ♦/*. 3. £ 




(2.5) 


because the stopping power is 


¥ £ > - !< £ > e , 

n 


( 2 . 6 ) 


and the atomic cross section is 


oj> l (£) = £<$(£) ( 2 - 7 ) 

n 

Equations (2.5)-(2.7) permit the rewriting of equation (2.3) in the usual continuous slowing-down 
approximation as 

£2 • V<|>.(i, £2, E)-l[Sj{E)*jl, £2, eJ] +a / (E)^(i. £2, eJ 

= J ’ £ ’ d£ (2,8) 

k 

where the cross section of equation (2.8) now contains only the nuclear contributions. 
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3. Transport Formalism 

The Boltzmann equation (2.8), as derived in section 2, can be rewritten as 

*' a e) = £J 'dE da a Jk [E, E, a, Q , e) (3.1) 

k 

where tyj^x, Q, eJ is the flux of ions of type; with atomic mass A. at x with motion along Q and 
energy E in units of A MeV, Cj(E) is the corresponding macroscopic cross section, 5\ ( E) is the linear 

energy transfer (LET), and cr^E, £ , Q, a J is the production cross section for type j particles with 

energy E and direction a by collision of a type k particle of energy E and direction ft'. The term that 
contains Sj (£) on the left side of equation (3.1) is the result of the continuous slowing-down approxi- 
mation. The solutions of equation (3.1) are unique in any convex region for which the inbound flux of 
each particle type is specified everywhere on the boundary surface. If the boundary is given as the loci 

of a two-parameter vector function y (s, t) for which a generic point on the boundary is given by f, 
then the boundary condition is specified by requiring the solution of equation (3.1) to meet 

%{ r , a, e) = f, a, e J (3.2) 

for each value of Q such that 

Qn[rj<0 (3.3) 

where n[rj is the outwardly directed unit normal vector at the boundary surface at point f, and \|/ is 
the specified boundary condition. 

The fragmentation of the projectile and target nuclei is represented by the quantity 
o^Zs, E , Q, Q, J , which is composed of three functions 

° jk ( E • a O J = O k (£') m jk {E)f.^E, £', Q, Q'J (3.4) 

where m {E) is the multiplicity (i.e., average number) of type j particles being produced by a colli- 
sion of type k of energy E , and E , Q, Q J is the probability density distribution for producing 

particles of type j of energy E in the direction £1 from the collision of a type k particle with energy E 

moving in the direction Q . For an unpolarized source of projectiles and targets, the energy angle distri- 
bution of reaction products can be a function of energies and cosine of the production angle relative to 
the incident projectile direction. The secondary multiplicities ( E ) and secondary energy angle dis- 
tributions are the major unknowns in ion transport theory. 

Information on the multiplicity m. k ( E ) was obtained in the past through experiments with galactic 

cosmic rays (GCR) as an ion source, and the fragmentation of the ions on target nuclei was observed in 
nuclear emulsion. (See ref. 23.) Such data are mainly limited by not knowing precisely the identity of 
the initial or secondary ions and by relatively low-counting rates of each ion type. The heavy ion accel- 
eration by machine reduces the uncertainty because high-counting rates can be obtained with known ion 


'^-x.TE s A E >+ a A E > 


V 
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types. In addition, accelerator experiments provide information, which was not previously available, for 
the spectral distribution E, ft, ft J. (See ref. 24.) 

The spectral distribution function consists of two terms that describe the fragmentation of the 
projectile and the fragmentation of the struck nucleus as follows (refs. 25 and 26): 

O jk [E, E , £2, Q'J = c k (E) [vj (E')/f k [E, E, £2, ft J + vj k (E^fJ^E, E, Q, fl'J] (3.5) 

where v£ and f? k depend only weakly on the target and vj k and fj k depend only weakly on the projec- 
tile. Although the average secondary velocities associated with f* are nearly equal to the projectile 
velocity, the average velocities associated with f T are near zero. Experimentally, 



m 

3/2 

JlE exp 

1 


. 2ic (°/*) 2 




m 

3/2 

JlE exp 

( aj2mE-aj2mE ‘ j 


2n (° 5) 2 

21 



(3.6) 


where p and p are the momenta per unit mass of j and k ions, respectively, and 

1 3/2 

J2E exp | 


fj k { E, E'.Q.Q' J- 


m 








(3.7) 


where aj k and aj k are related to the rms momentum spread of secondary products. These parameters 
depend only on the fragmenting nucleus. Feshbach and Huang (ref. 27) suggested that the parameters 
c? k and aj k depend on the average square momentum of the nuclear fragments as described by Fermi 
motion. A precise formulation of these ideas in terms of a statistical model was obtained by Goldhaber. 
(See ref. 28.) 

4. Approximation Procedures 

4.1. Neglect of Target Fragmentation 

The use of equations (3.5)-(3.7) in the evaluation of the source term t, 1 i, £2, £ J of equation (3.1) 
results in 

£,(*, a, e) = ZjdE' da G k (£') * k [i, £2 , E’ ) [v£ {E)fr k {E, E, £2, £2'J + vj fJ^E, E, £ 2 , fij] 

k 

= ^(U,£) + ^U,£) (4 ‘ 1} 

where as before, the superscripts P and T refer to the fragmentation of the projectile and the target, 
respectively. The target term is 


y(l, £2, e) = £ 


m 


2n\ aj k 


3/2 


J2E 


exp 


mE 


w 


jdaf dE vj (E) G k (E) O, e) (4.2) 
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which is negligibly small for 



m 

Thus, for calculating the flux at high energy, 

a n, e) 


(4.3) 


(4.4) 


4.2. Space Radiations 

Space radiations have the convenient property of being nearly isotropic. This fact, coupled with the 
forward-peaked spectral distribution, leads to a substantial simplification in the source term as follows: 


e) - 'ZjdE' dn G k (£') v£ (£') 


2rt(CTp' 


3/2 


4ie' 


exp 


(flV2>n£-flV2m£') 




t (in',£') 


(4.5) 


■4 / 


If £2 , £ J is assumed to be a slowly varying function of £2', an expansion about the sharply 
peaked maximum of the exponential function is possible. Such an expansion is made by letting 


£2 = Q+ ( cosG - 1 ) £2 4 - e^sinB 


(4.6) 


where 


and 


The flux may be expanded as 


cos0 = Q £1 


f 

> £2 x £2 

e <l> “ -» 

£2 x £2 


(4.7) 


(4.8) 


i(*. “'• / - 


E I + 




- [ (cosG - + sin0] + ... 


Substitution of equation (4.9) into equation (4.5) and simplification result in 


(4.9) 


k 


m 


2n 


P y 

J k J 


3/2 


A 

Je' 


x j <t> .( X, £2, E - 


£2 • A A l a e) 

d£2 V ' 


exp 


Jk 


C JlmE - JlmE 


2mjEE 


(4.10) 
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The leading term of equation (4.10) is clearly a good approximation of the source term whenever 


2mE 




(4.11) 


The leading term is equivalent to assuming that secondary ions are produced only in the direction of 
motion of the primary ions. For space radiations, which are nearly isotropic, equation (4.1 1) is easily 
met, and neglection of higher order terms in equation (4.10) results in the usual straight-ahead approxi- 
mation. If radiations are highly anisotropic, then equation (4.11) is not likely to apply. Validity of 
straight-ahead approximation was studied empirically by Alsmiller et al. (refs. 29 and 30) for proton 
transport. 


4.3. Velocity-Conserving Interactions 

Customarily, in cosmic ion transport studies (ref. 31), fragment velocities are assumed to be equal 
to the fragmenting ion velocity before collision. The order of approximation resulting from such an 
assumption is derived with the assumption that the projectile energy E is equal to the secondary energy 
plus a positive quantity e, 

E' = E + e (4.12) 

where e is assumed to contribute to equation (4.10) only across a small range above zero energy. 
Substitution of equation (4.12) into equation (4.10) and expansion of the integrand result in 





(E)vf k (E) 


F d Ai 

£ dE +* 


X , £ 2 , 



Q, E 


a 


jk 


2 mE 


(4.13) 


Because /mE« 1, the assumption of velocity conservation at those energies for which most 

nuclear reactions occur is inferior to the straight-ahead approximation but may be adequate for space 
radiations where variation of <!>*(-*, £2. £ J with energy is sufficiently smooth. That is 





£2, E 


£ 2 , £■] 


4.4. Decoupling Target and Projectile Flux 

Equation (3.1) can be rewritten with equation (4.1) as 

Bj$j{x, £2, e) = Z F J kh [l £2, e) + O* E ) 

k k 

where the differential operator is given by 

3 



(4.14) 


(4.15) 
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and the integral operator F. k = Fj k + F? k is given by 


F jkh{x> e) = JdE dQ c.^E, E, ft, ft )<t>^, ft , e) 

(4.16) 

By defining the flux as the sum of two terms, 


4 >/(•*> £2, e) + ft, e) 

(4.17) 

which permits the following separation: 


*/>/[*• “• e ) = 

k k 

(4.18) 

b j*A j - “• E ) = 4. eJ + 5, eJ 

k k 

(4.19) 

As noted in connection with equations (4.1)— (4.4), the source term on the right side of equation (4.19) is 
small at high energies. Assume that 

<!>/(*, ft, eJ * o 

(4.20) 

for E » ^CJj^ j / m . As a result of equation (4.20) and the fact that the ion range is small compared with 
its mean free path at low energy, 

k 

(4.21) 

E ) * E ) 

(4.22) 


k 


The advantage of this separation is that, once equation (4.21) is solved, equation (4.22) can be solved in 
closed form. The solution of equation (4.22) is accomplished by noting that the inwardly directed flux 
<| must vanish on the boundary so that 


<t> 


]{x, 3, E 


-zf 


dE 


A .P (E) 

J J_ 

Pj(E)Sj(E). 


jdE" dn ojjyE, e", ft, ft') 


X {x+ [Rj(E) -E.(E)]ft, ft, E"} (4.23) 

where = Rj^ [d + R. ( E ) ] with d being the projected distance to the boundary. 

The use of equations (3.5) and (3.7) in equation (4.23) yields 


where 


*/(U. e ).^ 


dE 


A.P.(E') 

J J K ’ 

m 

Pj(E)S.(E) 



13/2 


Jie' 


exp 


mE 

Wf 


x;/{x+ [E.(E) -E.(£')]ft} 


(4.24) 


5/ (*) = X \ dE a k ( £ ) v Jk ( £ ) O'. E ') (4.25) 

k 
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and o\ has been assumed to be a slowly varying function of projectile type k and projectile energy E. If 
the range of secondary type j ions is small compared with their mean-free-path lengths and the mean 
free paths of the fragmenting parent ions l k , then 


2*1 


m 


« l,. 


(4.26) 


and the integral of equation (4.24) may be simplified as 


tK** a ' ~ Sj (E) £ 


m 


3/2 


/ \2 aJ'ZE CXp / \o 

L (°s)1 


mE 


dE 


(4.27) 


which can be reduced into terms of known functions. Thus, 

A. 




E I = 


JrT ( > 


Sj(E) 


qCx) 


1 


2njn 


3 mE 
2 ( t Y ~ r 


3 mE y 


HI 


in terms of the incomplete gamma function. Equation (4.28) can be shown to be equivalent to 

A . 




+ 


mE 


71 a 


jk 


exp 


mE m Ey mEo 


At points sufficiently far from the boundary such that 

RJ X (d) » \ j ' 


m 


equation (4.28) may be reduced to 


ft, 


Hi 


Sj ( E) V ( 2n}2 


erfc 


mE I mE mE 


(4.28) 


(4.29) 


(4.30) 


(4.31) 


which is the equilibrium solution because the target fragment spectrum is the difference between the 
collision source and collision losses. The solution of equation (4.21) is examined in section 4.5. 


4.5. Back Substitution and Perturbation Theory 

One approach to the solution of equation (4.21) results from the fact that the multiple charged ions 
tend to be destroyed in nuclear reactions. Thus, 

Fj k = 0 0 > *) < 4 - 32 ) 

This means that there is a maximum j such that 

ft, e) = 0 (4.33) 
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where J is the maximum j. Furthermore, 

- fT-l.XUS-fJ (4.34) 

and in general. 


B 


'-tftj- ^ h F J-N,J - 




,p fi 

J-k \ x ’ 


ft, E 


(N<J- 1) 


(4.35) 


*= l 

Note that equations (4.34) and (4.35) constitute solvable problems. The singly and doubly charged ions 
satisfy 


~ F \, £ ’)+ ^ f] 


(4.36) 


k = 2 


Equation (4.36), unlike equations (4.33)-(4.35), is an integral-differential equation that is difficult to 
solve directly. Equation (4.35) is solvable by perturbation theory, and the resultant series is known to 
converge rapidly for intermediate and low energies. (See refs. 10 and 32-34.) Note that equations (4.33) 
and (4.35) are also obtained from perturbation theory as applied to equation (4.21) at the outset. Thus, 
the perturbation series is expected to converge after the first J plus a few terms. 


5. Galactic Ion Transport 


5.1. Derivation of GCR Heavy Ion Transport Equation 

In this section, the methods of previous nucleon transport studies (ref. 32) are expanded by combin- 
ing analytical and numerical tools. The galactic cosmic ray ion transport problem is transformed to an 
integral along the characteristic curve of that particular ion, and the perturbation series (ref. 32) is 
replaced by a simple numerical procedure. The resulting method reduces the difficulty associated with 
the low-energy discretization and the restriction to a definite form for the stopping power. The resulting 
numerical procedure is simple and is not computationally intensive. 

Here, the straight-ahead approximation is used, and the target secondary fragments are neglected. 
(See refs. 6 and 8.) For multiple charged ions, the transport equation can be written as 


ra 

.dx 


Ye s J (E)+ °j] { x ’ E) = X m ,*° A 

k>j 


(5.1) 


where <|>y (x, E) is the flux of ions of type j with atomic mass A . at x moving along the X-axis at 
energy E in units of A MeV , cr is the corresponding macroscopic nuclear absorption cross section, 
(E) is the change in E per unit distance, and nij ^ is the multiplicity of ion j produced by ion A: in a 
collision. The corresponding transport equation for the light ions is 


rd_ 

_9x 


+ V £) ]v*’ E) = If °j k (^E')<\> k (x,E') dE 


(5.2) 


The quantities m jk and o y are assumed to be energy independent in equation (5.1) but are fully energy 
dependent in equation (5.2). 


The range of the ion is given as 



(5.3) 
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The solution to equation (5.1) is found subject to the boundary specification at x = 0 and arbitrary E as 


♦ y (0 ,E) = F.{E) 


Usually, Fj(E) is called the incident beam spectrum. 
From Bethe’s theory. 


Sj(E) 


A 7 2 

a£’ ie> 

J p 


(5.4) 


(5.5) 


which holds for all energies greater than =100A keV, provided that the ion effective charge is used, and 
leads to 


Z 2 Z 2 

J-Rj(E) = fR p (E) 

J P 


(5.6) 


The subscript p refers to proton. Equation (5.6) is accurate at high energy and only approximately true at 
low energy because of electron capture by the ion (which effectively reduces its charge), higher order 
Bom corrections to Bethe’s theory, and nuclear stopping at the lowest energies. Herein, the 
parameter Vj is defined as 

(5.7) 


V.R.(E) = v,* 4 (£) 


so that 


Z? 


v = -r- 

J A. 

J 


(5.8) 


Equations (5.6M5.8) are used in the subsequent development, and the energy variation in Vj is 
neglected. The inverse function of Rj ( E ) is defined as 


E = RjHRj(E)] 

For the purpose of solving equation (5.1), the following coordinates are used: 

ri jmx-RjiE) 

+ Rj (E) 

where X] . varies along the particle path, and £.. is constant along the particle trajectory. The new fluence 
functions are defined as 


(5.9) 

(5.10) 

(5.11) 


where 


Xj (Oy. kj) =Sj(E) 4> y (*, E) = \\ h (*, r.) 

S k ) 

*V + T 1 j = ^ + T li 


’irt = 7 ( V5t) 

j 


(5.12) 

(5.13) 

(5.14) 

(5.15) 
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(5.16) 


r. - R. (E). By this coordinate mapping, equation (5.1) becomes 

( 2 4' + 0 /) £ ; < ’VV = 

J /. k 

where Oy is assumed to be energy independent for simplicity of formalism. There is a small variation 
in o. (=20 percent) that is accounted for in the code. (See ref. 35.) Equation (5.16) is solved by using 
line integration with an integrating factor (ref. 12), which results in 



*v (T vV = ex p[k^ + v] 

(5.17) 

Then, 

Z y (tl /V = eX p[-^(4y + Tly)]XyK / ^) 




(5.18) 

where 

, V * + V / - V i“ V / 



V *- V ) , V * +V /t 
- 2v t " + 2v, % 


By defining 

o 

II 

s* 

s? 

(5.19) 

then 




Vj (*> r ;) _ ex P Vy ( 0. r j + x) + J * dz exp {-ajz) ^ m jk 0 k -^W k x~z,r k + -+z 

0 * V * V v k . 

Furthermore, 

V|T (x + h, rj) = exp (-oTi) i ]f. (x, r. + h) 

r ' 1 v. f v. > 

+ dz exp (-oz) X m jk c k^k x + h -t’ r k + Z L z 

u k k V v k J 

Equation (5.20) clearly shows that 

\\t k (x + h-z,r k ) = expl-a k (h-z)]y k (x,r k + h) +0{h-z) 
which, after substitution into equation (5.21), yields 
Vj ( x + h, r.) = exp ( -G.h ) \|r (x, r. + h) 

di _ V . v . 

+ \dz exp (-OjZ) X m jk a k^f ex P [-O* (A -z) ] V* X, r* + -l z + h - z 
0 *• * v * 


(5.20) 


(5.21) 


(5.22) 


(5.23) 
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which is correct to order h 2 . This expression can be further approximated by 


\| f.(x + h,rj) = exp(-Gjh)y^j(x,rj + h) 


v . I" exp (-a.A) - exp (-O^A) 
a k~°j 


1 1 

V. > 

> 

1 

r k + ^ h 

\ v * J 


(5.24) 


which is accurate to O [ (v t - v .) A]. Equation (5.24) is the basis of the GCR transport computer code 
GCRTRN. (See refs. 12, 13, andf 36.) The nucleon transport equation (5.2) is solved by adding the heavy 
ion collision source of nucleons to the BRYNTRN computer code (ref. 33) to create HZETRN, which 


effectively solves equation (5.2) by adding a heavy ion collision source of nucleons to the right side of 
the equation. Equation (5.24) provides the propagation algorithm for the heavy ions. The corresponding 
propagation procedure for the light ions (refs. 33, 37, and 38) is 


\\>(x + h,r) ~ e~ ah \\i (x, r + h) + e-° h f dzf dr'f(r + z, r' + z)\ V (*, r + h) 

JQ J r 


(5.25) 


with the order of h 2 . 

The following quantities are of interest: 

1 . The integral fluence is 

V*- >£) - * 


(5.26) 


2. The energy absorbed per gram is 


D (x,>E) = f° A. \|/.[x, RAE)] dE 

J JR J J 3 

3. The dose equivalent is 


(5.27) 


H. (x, > E) = Q F Vj [x, Rj (E) ] dE (5.28) 

where Q F is the quality factor. These quantities are not recommended for use in shield design for pro- 
tection from GCR but give some relative measure of their biological importance. 


5.2. Numerical Procedure 

The secondary particle production term of the propagation algorithm for the nucleons in 
equation (5.25) has been further reduced to a form that can be numerically integrated with ease. Details 
of the form and its validity are discussed in reference 33. For the heavy ions, the secondary production 
term (i.e., the second term on the right side of equation (5.24)) does not involve any integration; how- 
ever, the interpolation of the transformed fluence function is based on the independent variable r k , 
which is different from the range of ions r . of type j given on the left side of the equation. To circum- 
vent the problem, the equation is further modified by using the definition of Sj (E), with E given in units 
of A MeV; then 


Sj(E) = SjiE/Aj) 


Ax 


1 


1 AE i 

- J 1 = —S.(E) 

A. Ax A. J y r 
J J 


(5.29) 


with 


W = Z fV E J /A ? 


(5.30) 
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(5.31) 


where S p is the proton stopping power, and E. is the energy in MeV of the ions of type j. Then 


*/ <£) = z z j s p< e / a j> -V,<V 

J 


where 


E n = EVA . 

P J j 

By rewriting equation (5.19) as 

V, (-*> rj) = Sj (E) 4> 7 (X, E) = vS p (E) (x, E) 

the new fluence function can be defined as 


(5.32) 


(5.33) 


Vj (*> r ) - v j s p ( £ ) <t> (x, E) s y. (x, r.) 

with r = r p = \.r. , where 


Equation (5.24) becomes 


r 

P 


r 


d E 

o S p {E) 


(5.34) 


(5.35) 


Wj (x + h,r) = e 


- a.h 


J 'Vj(x,r + Vjh ) + £ 


m jk G k 




°7 


-\\f k (x,r + v.h) 


(5.36) 


so that there is now only one definition of range that is related to energy. The equation can be solved by 
setting up the proton range grid r and marching the solution from x = 0 by steps of h to the desired 
thickness. 


5.3. Error Propagation 

When considering how errors are propagated in the use of equation (5.36), the error is introduced 
locally by calculating \| f. (x> r + Vjh) across the energy range grid. By limiting the analysis to the first 
term of equation (5.36), the error is defined at each range grid r. such that 

Vj ( x + K r.) = e~ a j\j (x, r. + v.h) (5.37) 

The truncation error e. is introduced in the interpolation procedure for the value \j / and 

Vj (x, r. + Vjh) = int (x, r. + Vjh) + e. (A) (5.38) 

After the wth step from the boundary, the numerical solution is 

m - 1 

r i) = ^VFj, int [ («- U*. ' , , + v / A] + £ e~°J im ~ l)h e,(h) (5.39) 

/ = 0 

Suppose 0 < (/z) < e (/i) for all /; then the propagated error is bounded by 

m - 1 m - l 

W" 11 = I sew I (5.40) 

I = o 1 = 0 
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Note that 


m - 1 

-O.mh O .hi 
2j e 1 e 1 

l = 0 




(5.41) 


because hOj « 1 . Clearly, the propagated error on the mth step is bounded by 


... e(h) 
e (h) < -p — - 
prp v ’ he. 


. -o.mh 

1 - e J 


(5.42) 


where e ( h ) is the maximum error per step. With the increasing value of m, the propagated error grows 
with each step to a maximum value of z{h) /hOj. Because the increase in the value of h is limited by 
the perturbation theory, the reduction of the local truncation error is the only viable approach left for 
reducing the propagated error to a desired level. The same consideration can be applied to the second 
term of equation (5.36) as the terms are similar in nature. The issue of error propagation in HZETRN is 
further studied by Shinn and Wilson. (See ref. 39.) 


5.4. Numerical Algorithms 

The error analysis of section 5.3 results in the conclusion that, to effectively reduce the level of 
propagated error, the local truncated error must be reduced. Three numerical algorithms are involved in 
solving equations (5.25)-(5.28) and (5.36): interpolation, integration, and grid generation. The choice of 
grid distribution that is interrelated to the interpolation and integration methods can increase the 
efficiency of the code if the number of grids can be reduced. 

The interpolation method in HZETRN is the third-order Lagrange’s method, which was used suc- 
cessfully in improving BRYNTRN. (See ref. 37.) With the four neighboring interpolation grids (data 
points) placed on both sides of the interpolated point, the error will tend to be the smallest in the middle 
interval of all the data points if the grid distribution is rather uniform. (See ref. 40.) The choice of a 
much higher order Lagrange’s method will substantially decrease the efficiency of the computer code 
because there are more than 10 interpolation calls for each energy point at each step. Other interpolation 
methods such as a cubic spline were considered but discarded. In general, the splines are more accurate; 
however, their characteristically large oscillations can result in erroneous solutions. 

The procedure for numerical integration that was used in the improved BRYNTRN (ref. 37) is also 
used here in HZETRN. The procedure is based on the compound quadrature formulation summing over 
all the subintervals between the grids with the midpoint evaluated by making use of the improved inter- 
polation procedure mentioned previously. A simple numerical method such as Simpson’s rule is used to 
integrate for the subintervals. 

Three considerations dictate how the grids should be distributed. The first is the shape of the input 
spectrum. Because the GCR fluences are several orders of magnitude greater at the lower energies 
(refs. 41—43), a logarithmic scale is used for the energy or range coordinate as was done in reference 37. 
The second is the choice of the interpolation method, which requires that the four neighboring grids be 
as uniformly spaced on a logarithmic scale as possible so that the interpolation error can be minimized. 
Because the interpolation is performed on the range grid rather than on the energy grid, a uniform grid 
distribution on a logarithm of range r is desired. The third is the code efficiency, which increased almost 
quadratically with the decrease in the number of grid points. With the uniform grid points on a logarith- 
mic scale r as the basic structure, the distribution can be modified further to reduce the number of points 
in the region in which the data are not propagating through the steps. For BRYNTRN, this ^region 
occurred below r mjn + h (i.e., =1 g/cm 2 ) because r mjn « 1 and h is assumed to equal 1 g/cm 2 . (See 
ref. 37.) This also^pplies to HZETRN, although the interpolation is now at ' min + v/i, where Vj is 
always equal to or much greater than 1 ; because in propagating a distance h, all ions with energy below 
r + h are slowed to energies below the first energy grid point. 

min 
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6. Stopping Power 

In passing through a material, an ion loses the greater fraction of its energy to electronic excitation 
of the material. Although a satisfactory theory of high-energy ion-electron interaction is available in the 
form of Bethe’s theory utilizing the Bom approximation, an equally satisfactory theory for low energies 
is not available. Bethe’s high-energy approximation of the energy loss per unit length (i.e., stopping 
power) is given by 


S 


e 


4nNZ 2 Z,e 4 

P 1 

2 

mv z 




( 6 . 1 ) 


where Z p is the projectile charge, N is the number of target molecules per unit volume, Z is the number 
of electrons per target molecule, m is the electron mass, v is the projectile velocity, {3 = v/c, c is the 
velocity of light, C is the velocity-dependent shell correction term (ref. 44), and 7 is the mean 
excitation energy given by 


Z / ln < 7 /) = X4 ln ( E „) (6.2) 

n 

wherc /„ are the electric dipole oscillator strengths of the target, and E n are the corresponding excita- 
tion energies. The sum in equation (6.2) includes discrete and continuum levels. Molecular stopping 
power is reasonably approximated by the sum of the corresponding empirically derived atomic stopping 
powers for which equations (6.1) and (6.2) imply that 

Z In (/) = £ nZ . In (/.) (6.3) 

j 

where Z and I pertain to the molecule, Z. and I. are the corresponding atomic values, and nj are the 
stoichiometric coefficients. This additive rule (eq. (6.3)) is usually called Bragg’s rule. (See ref. 45.) 

As an input to HZETRN, the high-energy cutoff for the incident GCR spectrum is usually taken to 
be far beyond the minimum stopping power of =24 GeV, where Bethe’s theory starts to overestimate 
and worsens as the energy increases. This overestimation can be corrected by considering that the field 
of the incoming ion projectile is perturbed by the neighboring polarized atoms. In HZETRN the approx- 
imate correction for the polarization effect (i.e., density effect) is based on the work of Stemheimer. 
(See ref. 46.) In reference 46, the polarization effect corrections for various elements and compounds 
were evaluated and fitted by the expressions 


8 = 0 

(x<x 0 ) 

(6.4a) 

8 = In r| 2 + C + a (xj -x) m 

(x g <x<x l ) 

(6.4b) 

8 = In q 2 + C 

(x>X|) 

(6.4c) 


where 8 is the density effect correction used in the stopping-power formula, q = p/m c , p is the 
momentum, m Q is the rest mass of the charged particle, c is the velocity of light, and 

x = (log 1Q e)ln r) = 0.43429 In q 

The quantities a , m, x q , and x, are constants that must be evaluated for each material; C is given by 

C = -2 In (I/hv p - 1) (6.5) 
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( 6 . 6 ) 


where / is the mean ionization potential; the plasma energy hv p of the material is 

hv p = /i|ne 2 /nm t j 

where n is the electron density in electrons/cm 3 , m g is the rest mass, and e is the charge of an electron. 
For compounds and mixtures, the effective mean ionization potential can be determined by 


In / = -Y n. In /. (6.7) 

i i 
t 

where n. is the electron density for the ith element, and /. is the corresponding atomic ionization 
potential. 

Reference 46 suggests that, for some practical applications, the use of only the asymptotic density 
effect correction (eq. (6.4c)) may be adequate for all charged particle energies. In reference 47 
Armstrong and Alsmiller compared the stopping-power differences between correct expressions and an 
asymptotic formula for several elements and compounds. The results indicated an overestimation of at 
most 6 percent when using the asymptotic formula. Therefore, only the asymptotic formula is used in 
HZETRN. Equation (6.1) with polarization effect modification becomes 
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4kNZ 2 Z e 4 
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2 mv 2 
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p2_C_8 } 

H Z t V 


( 6 . 8 ) 


The use of equation (6.8) in conjunction with the asymptotic formula of equation (6.4c) does not 
involve any evaluation of the constants a, m, x Q , and x, for the material considered. A comparison of 
equations (6.1) and (6.8) for different materials is discussed by Shinn et al. (See ref. 48.) 

The electronic stopping power for protons is adequately described by equation (6.8) for energies 
>500 keV for which the shell correction C makes an important contribution up to 10 MeV. (See ref. 49.) 
For proton energies <500 keV, charge exchange reactions alter the proton charge over much of its path 
so that equation (6.8) is understood to be an average over the proton charge states. Normally, an average 
over the charge states is introduced into equation (6.8) so that the effective charge is the rms ion charge 
and not the average ion charge. At any ion energy, charge equilibrium is established very quickly in all 
materials. By utilizing the effective charge in equation (6.8), only modest improvement results at ener- 
gies <500 keV, which presumably indicates the failure of this theory based on the Bom approximation. 
(See refs. 49 and 50.) 

The electronic stopping power for alpha particles requires terms in equation (6.8) of higher order in 
projectile charge Z because of corrections to the Bom approximation. The alpha stopping power can 
not be related to tlie proton stopping power through their effective charges. Parametric fits to experi- 
mental data are given by Ziegler in reference 51 for all elements in both gaseous and condensed phases. 

The electronic stopping powers for heavier ions are related to the alpha stopping power through 
their corresponding effective charges. HZETRN uses the effective charges suggested by Barkas 
(ref. 52) of 


Z* = Z p [l-exp[-125p/z2/3j] 
where Z p in equation (6.9) is the atomic number of the ion. 


(6.9) 
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At sufficiently low energies, the energy lost by an ion in a nuclear collision becomes important. The 
nuclear stopping-power theory used in HZETRN is a modification of the theory of Lindhard, Scharff, 
and Schiott. (See ref. 53.) The reduced energy is given as 

32.53 A A ( E 

£ = T 1 \wi ( 6 - w ) 

v> ( v^( z r +z < 2/3 ) 

where E is in units of A keV, and A p and A { are the atomic masses of the projectile and the target, 
respectively. The nuclear stopping power in reduced units (ref. 51) is 
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n 


1.59e 1/2 

1.7e 1/2 In (e + e 1 ) 
1 + 6.8e + 3.4e 3/2 
In (0.47e) 

2e 


(E<0.01) 

( 0.01 <£< 10 ) 

( 10 < e) 


( 6 . 11 ) 


and the conversion factor to units of ( e V-atoms/cm 2 )/l 0 1 5 is 

8.426Z A A 

f - P j P 

<v-v( z r +z ? /3 ) 1/2 

The total stopping power Sj is obtained by summing the electronic and nuclear contributions. Other 
processes of energy transfer such as Bremsstrahlung and pair production for stopping massive ions are 
unimportant. 

For energies greater than a few A MeV, Bethe’s equation is adequate if the appropriate corrections 
to Bragg’s rule (refs. 54-56), shell corrections (refs. 44, 49, and 50), and an effective charge are 
included. Electronic stopping power for protons is calculated from the parametric formulas of Andersen 
and Ziegler (ref. 49) with some modifications. (See ref. 9.) 

Because alpha stopping power is not derivable from the formula for proton stopping power by using 
the effective charge at low energy, the parametric fits to empirical alpha stopping powers given by 
Ziegler (ref. 51 ) are used. Physical state and molecular binding effects are most important for hydrogen 
(ref. 54), water stopping power was approximated by using the condensed phase parameters for hydro- 
gen and the gas phase experimentally derived parameters for oxygen. Electronic stopping powers for 
ions with a charge greater than 2 are related to the alpha stopping power by the effective charge of equa- 
tion (6.9). Figures 2(a), 2(b), and 2(c) show the transport coefficient ranges and stopping powers for five 
different charge values Z, as calculated by HZETRN for aluminum, water, and liquid hydrogen targets, 
respectively. 


7. Nuclear Database 

The nuclear cross sections for neutron and proton interactions are described extensively in refer- 
ences 10 and 33. The heavy ion absorption cross sections o abs are currently derived from 

°abs = Kr o + Aj/ 3 - ^0.200 + y4~' - 0.292 exp^- jcos (0.229£°- 453 ) J } 2 ( 1 + 5£~') (7.1) 

which for r Q = 1.26 fm were fit to the quantum mechanical nuclear cross sections calculated in 
reference 57. 

In this section, the heavy ion fragmentation mechanism is described by abrasion-ablation formal- 
ism. In the abrasion-ablation fragmentation model, the high-energy projectile nuclei collide with 
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stationary target nuclei. In the abrasion step, those portions of the nuclear volumes that overlap are 
sheared away by the collision. The remaining projectile piece, called “a prefragment,” continues its tra- 
jectory with essentially its precollision velocity. As a result of the dynamics of the abrasion process, the 
prefragment is highly excited and subsequently decays by the emission of nuclear particles. This step is 
the ablation stage. Final deexcitation is through gamma emission. The resultant isotope is the nuclear 
fragment whose cross section is measured. The abrasion process can be analyzed with classical geomet- 
ric arguments (refs. 11, 36, and 58) or methods obtained from formal quantum scattering theory. (See 
refs. 59 and 60.) The ablation stage can be analyzed with geometric arguments (ref. 1 1) or more sophis- 
ticated methods based on Monte Carlo or intranuclear cascade techniques. (See refs. 61-64.) Predictions 
of fragmentation cross sections on hydrogen targets are made with the approximate semiempirical 
parameterization formulas of references 65 and 66. The fragmentation cross sections for other heavier 
targets are generated by the NUCFRG series of semiempirical fragmentation codes, in particular, the 
second version NUCFRG2 as described in reference 67. 

The amount of nuclear material stripped away in the collision of two nuclei of radius R p and R T is 
taken as the volume of the overlap region times an average attenuation factor. (See refs. 36 and 67.) The 
relevant formula for the constituents in the overlap volume in the projectile is given by 


A abr 

= FA p [ 1 — exp (-Cp/X) ] 

(7.2) 

where C T is the chord length at maximum overlap density of the intersecting volume of the projectile 
and the target X is the nuclear mean free path, and the expression for F depends on the nature of the col- 
lision (i.e., peripheral versus central) and the relative sizes of the colliding nuclei. 

From reference 61, for R T > R p , 



P = 0.125 (pv) 1/2 (i-2)( 

^) 2 -0.125[0.5(pv) 1/2 (i-2) + l] 

( 7 ) (7 - 3) 

and 



F = 0.75(1 -v) 1/2 ( 

0.125 [3(1 v) 1/2 ~ 1] ( 

f (7.4) 

with 

R P 

v = 

R p 4- Rj> 

(7.5) 


p= * 

Rp + Rj, 

(7.6) 

and 

1 . r t 

^v- 1 =R f 

(7.7) 


Equations (7.3) and (7.4) are valid when the collision is peripheral (the two nuclear volumes do not 
completely overlap). In this case, the impact parameter b is restricted such that 

R T -R p <b<R T +R p (7.8) 
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If the collision is central, then the projectile nucleus volume completely overlaps the target nucleus vol- 
ume ( b < R T -R p ), and all the projectile nucleons are abraded. In this case, equations (7.3) and (7.4) are 
replaced by 


and 


F = 1 

and there is no ablation of the projectile because it was destroyed by the abrasion. 

For a peripheral collision with R p >R T , equations (7.3) and (7.4) (ref. 63) become 

P = 0.125 (HV) 1/2 (‘ j 2 -Q.125{0.5^] 1/2 ^l 2 ^ (1/v) ^ ~M 2 ) 1/2 - H [ (2 - ^i] 1/2 ^ 1 - g J 

and 

F = n7sn-v^ H - d -p 2 ) 3/2 l [1 - (1 -u) 2 l l/2 ; |^1 -3j 3 

where the impact parameter is restricted such that 


(7.10) 


(7.11) 


(7.12) 


R p -R T <b<R p + R T 


(7.13) 


For a central collision (b<R p -R T ) with R p >R T , equations (7.1 1) and (7.12) become 


and 





1/2 


F = 


[1- (1 ~H 2 ) 3/2 J 



1/2 


(7.14) 


(7.15) 


The charge ratio of the nuclear matter removed by the collision is assumed to be that of the parent 
nucleus. The spectators in the overlap region are assumed loosely bound to the remaining prefragment 
and are lost in a preequilibrium emission process. (See ref. 67.) 

The surface distortion excitation energy of the projectile prefragment following abrasion of the 
m nucleons is calculated from the clean-cut abrasion formalism of reference 58. For this model, the col- 
liding nuclei are assumed to be uniform spheres with radii R. (i = P, T). In collision, the overlapping 
volumes shear off so that the resultant projectile prefragment is a sphere with a cylindrical hole gouged 
out of it. The excitation energy is then determined by calculating the difference in surface area between 
the distorted sphere and a perfect sphere of equal volume. This excess surface area AS is given in 
reference 61 as 


A5 = 4tc/? 2 [ i +P _ (i _F)2/3] (7.16) 

where the expressions for P and F differ and depend upon the nature of the collision (i.e., peripheral 
versus central) and the relative sizes of the colliding nuclei, which were given in previous equations. 

The excitation energy associated with surface energy is known to be 0.95 MeV/fm 2 for near- 
equilibrium nuclei so that 


E' s = 0.95AS 


(7.17) 
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for small surface distortions. When large numbers of nucleons are removed in the abrasion process, 
equation (7.17) is expected to be an underestimate of the actual excitation. This requires the introduc- 
tion of an excess excitation factor (refs. 1 1 and 36) in terms of the number of abraded nucleons A abr as 

/ = 1 +^L r + [1500-320(A- 12)]^ (7.18) 

p A * 

P P 

which approaches 1 when the impact parameter is large but increases the excess excitation when many 
nucleons are removed and grossly distorted nuclei are formed by the collisions. (See ref. 11.) The quan- 
tity in the brackets accounts for light nuclei having greater excitation energy than heavier nuclei for the 
same fraction of mass removed (ref. 67) and is limited to values between 0 and 1500. The total excita- 
tion energy is then 


E = 

s 


V 


(7.19) 


which reduces to equation (7.17) for small A abr . Also, the assumption is that 100 percent of the frag- 
ments with a mass of 5 are unbound, 90 percent of the fragments with a mass of 8 are unbound, and 
50 percent of the fragments with a mass of 9 ( 9 B) are unbound. (See ref. 11.) 

A secondary contribution to the excitation energy is the transfer of kinetic energy of relative motion 
across the intersecting boundary of the two ions. The rate of energy loss of a nucleon when it passes 
through nuclear matter (ref. 68) is taken as 13 MeV/fm, and the energy deposit is assumed to be sym- 
metrically dispersed about the azimuth so that 6.5 (MeV/fm)/nucleon at the interface is the average rate 
of energy transfer into excitation energy. (See ref. 36.) This energy is transferred in single-particle colli- 
sion processes, and the energy is transferred to excitation energy of the projectile for half of the events 
and leaves the projectile excitation energy unchanged for the remaining half of the events. (See ref. 1 1 .) 
The first estimate of this contribution uses the length of the longest chord Cj in the projectile surface 
interface. This chord length is the maximum distance traveled by any target constituent through the pro- 
jectile interior. The number of other target constituents in the interior region can be found by estimating 
the maximum chord C { transverse to the projectile velocity, which spans the projectile surface inter- 
face. (See ref. 36.) The total excitation energy from excess surface and spectator interaction is then 


E = 13 C, 

X 1 


l+i(C,-1.5) 


(7.20) 


where the second term contributes only if C f > 1 .5 fm . Also, the effective longitudinal chord length for 
these remaining nucleons is assumed to be one third of the maximum chord length. 

The decay of highly excited nuclear states is dominated by heavy particle emission. In the present 
model, a nucleon is assumed to be removed for every 10 MeV of excitation energy. The number of 
ablated nucleons is 


A abl “ 


E S + E X 
10 


(7.21) 


In accordance with the previously discussed directionality of the energy transfer, E % has two values 


E = 1 

X ' 


E for P = ^ 
x x 2 


0 for E x = 2 


(7.22) 


where P x and P. are the corresponding probabilities of occurrence of each value in collisions. 
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The number A A of nucleons removed through the abrasion-ablation process is given as a function of 
the impact parameter as 


M =4 tbr (4) + i ib| (i,) 


The cross section for the removal of A A nucleons is estimated as 


(7.23) 


O(AA) = n\^b*-b^ 


(7.24) 


where is the impact parameter for which the volume of the intersection of the projectile contains 
^abr nuc l eons an d the resulting excitation energies release an additional A ^ nucleons at the rate of 
1 nucleon for every 10 MeV of excitation such that a ' 


similarly, for b j 


^abr^l) + ^abl^2^ “ ^ “ 2 
^abr + ^abl ~^ + 2 


(7.25) 


(7.26) 


In the previous discussion, the assumption of straight-line trajectories makes the impact parameter b 
the distance of closest approach. This assumption makes the abrasion-ablation model a high-energy 
heavy ion fragmentation model; thus, Coulomb corrections due to low-energy trajectories must be 
added to the standard abrasion-ablation model. Equations (7.24)-(7.26) can be combined and rewritten 
(ref. 58) as 


ct(AA) 



(7.27) 


Equation (7.27) is retained in NUCFRG2 as written. However, b is no longer taken to be the dis- 
tance of closest approach as used in equation (7.27) but is redefined in the following discussion. 

The equations of motion in the nuclear Coulomb field are given by energy conservation as 


E tot = + 


ZpZji ? 2 


2 fir 2 


(7.28) 


where £ tol is the total kinetic energy in the center-of-mass system at great relative distances, r is the rel- 
ative distance between the charge centers with time derivative r , p is the reduced mass, / is the angular 
momentum, Zp and Zj are the atomic numbers of the projectile and target nucleus, respectively, and e 
is the electronic charge. The angular momentum is given as 


/2 = 2 ^ E toi b2 (7.29) 

With the use of equation (7.29) and r = 0, equation (7.28) becomes 


which can be written as 


P _ *tot* 2 ( Z/> Z 7^ 
rot 9 + 


b 2 - r(r - r ) 

v nr 


(7.30) 


(7.31) 
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where 


r 


m 


ZpZj£^ 


(7.32) 


Note that r is the distance of closest approach for the zero impact parameter. For a given A A, 
equation (7.32) is used in NUCFRG2 to calculate the distance of closest approach and equation (7.31) is 
used to calculate the impact parameter, which permits extrapolation backward in time along the 
Coulomb trajectory to the initial impact parameter. (See ref. 69.) This calculated value of b is used in 
equation (7.27) to evaluate the cross section. Note that the effect of the Coulomb trajectory is to move 
the separation at impact r to smaller impact parameters b and thus reduce the cross sections, especially 
at low energy. 

A second correction to the trajectory calculation comes from the transfer of kinetic energy into 
binding energy in the release of particles from the projectile. The total kinetic energy in passing through 
the reaction zone is reduced from the initial energy E . to the final energy Ej 

E f =E.-\OAA (7.33) 

by assuming that 10 MeV is the average binding energy. The kinetic energy used in the closest approach 
calculation is the average 


£ tot = \( £ i + E ? = V^IOAA) < 7 - 34 ) 

As given by equation (7.34), E is obviously very crude and equation (7.34) can be refined further. 

The distribution of charge Z p of final projectile fragments of mass A p are strongly affected by 
nuclear stability. The charge distribution for a given a(AA) as described by Rudstam (ref. 70) is 

aiA^Zp.) = F^expy-R^Zp-SAp+TAp^ Ja (AA) (7.35) 

where R = 1 1 ,8/A j? , D = 0.45 , S = 0.486 , T = 3.8 x 10" 4 , and F, is a normalizing factor such that 

^c(A r Z p ) = o(AA) (7.36) 

The Rudstam formula for o(AA) was not used because the A A dependence is too simple and breaks 
down for heavy targets (refs. 11, 36, 71, and 72) and was replaced by equation (7.27) in the 
abrasion-ablation fragmentation model. 

The charge of the removed nucleons A Z is calculated according to charge conservation 

Z p = Z p + AZ (7.37) 


and is divided between the nucleons and hydrogen and helium isotopes according to the following four 
rules: 

1. The abraded nucleons are those removed from that portion of the projectile in the region of over- 
lap with the target. (See ref. 67.) Therefore, the abraded nucleon charge is assumed to be proportional to 
the charge fraction of the projectile nucleus as 


Z 


abr 


Z P A &br /A P 


(7.38) 
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This, of course, ignores the charge separation caused by the giant dipole resonance model of 
reference 63. The charge release in the ablation stage is then 

Z abl = AZ - Z abr (7-39) 

which simply conserves the remaining charge. 

2. The alpha particle is known to be unusually tightly bound in comparison with other nucleon 
arrangements. Because of this unusually tight binding of the alpha particles, the helium production is 
maximized in the ablation process. (See refs. 1 1 and 36.) The number of alpha particles is 

= f Int ( Z abl /2 ) ’ Int (^ab/ 4 ) 1 min (7.40) 

where Int(;t) denotes the integer part of x. The remaining isotopes are likewise maximized from the 
remaining ablated mass and charge in the order of decreasing binding energy. (See ref. 69.) The number 
of protons N produced is given by charge conservation as 

N P = Z abl-S Z ^, (7.41) 

i 

where i ranges over all ablated particles of mass 2, 3, and 4. Similarly, mass conservation requires the 
number of neutrons N produced to be 


N n = A A\- N p-'L A i N i < 7 - 42 ) 

i 

3. The calculation is performed for A A — 1 to A A — A p — 1, for which the cross section associated 
with A A >A^-0.5 is omitted. These are, of course, the central collisions for which the projectile is 
assumed to disintegrate into single nucleons if r p < r T . Then 

N P = Z /> (7.43) 


N = 

n 


(7.44) 


Other collision products such as the energetic target fragments and mesonic components are ignored. 
The peripheral collisions with A A < 0.5 are also ignored. Most important in these near collisions is the 
Coulomb dissociation process. (See ref. 73.) 

4. The description of the nuclear radius is given by 


R =' 29 l R L>-K ) (7-45) 

where is the rms radius of the nucleus charge distribution and R is the radius of the proton charge 
distribution. The improvements to the abrasion-ablation model as implemented in NUCFRG2 and com- 
parison with experiments are further described in references 74 and 75. 

Finally, the electromagnetic (EM) dissociation cross-section contributions to the nuclear database 
as implemented in NUCFRG2 must be considered. In electromagnetic dissociation, the virtual photon 
field of the target nucleus interacts electromagnetically with the constituents of the projectile to cause 
excitation and eventual breakup. The electromagnetic dissociation model implemented in NUCFRG2 is 
limited to single nucleon (proton or neutron) removal processes. 
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The total electromagnetic cross section for one nucleon removal resulting from electric dipole (£1) 
and electric quadrupole (£2) interactions is written 


°em ~ °E 1 + °E2 

= J[N £ 1 (£)o £1 (£) +N E2 (E)c e 2 (E)] dE 


(7.46) 


where the virtual photon spectra of energy £ produced by the target nucleus (ref. 76) are given for the 
dipole field by 


"bi 


= \-Z 2 a\ 

En |3 2 




(7.47) 


and for the quadrupole field by 


N e2 (E) 


1 2 2 1 

= --Z 1 a— 
(5 4 


En 


2 ( 1 - p 2 ) K\ + \ (2 - p 2 ) 2 K q K x - \w{k 2 - K 2 


(7.48) 


The terms a £1 (£) and a E2 (E) are the corresponding photonuclear reaction cross sections for the 
fragmenting projectile nuclei. The terms K {) and K x in the expression for N £] and N E2 are modified 
Bessel functions of the second kind and are also functions of the parameter %. The latter is 


\ = 


InEb. 

mi 

yP/ic 


(7.49) 


where E is the virtual photon energy, b min is the minimum impact parameter below which the collision 
dynamics are dominated by nuclear interactions rather than by EM interactions, p is the speed of the tar- 
get measured from the projectile rest frame as a fraction of the speed of light c, h is Planck’s constant, 
and y is the usual Lorentz’s factor from special relativity y - (1 -P 2 )” 1/2 . The minimum impact 
parameter is 


b . = ( 1 + * ,) b + — - 

min v a' c 271 


where x , = 0.25 and 


ZpZj £ 2 


w o v 


(7.50) 


(7.51) 


allows for trajectory deviation from a straight line. (See ref. 77.) The critical impact parameter for 
single-nucleon removal is 

b c = 1.34 [_Ap /2 +/4 j/ 3 - 0.75^Ap 1/3 +/f £ 1 / 3 J] (7.52) 

with A p and A T being the projectile and target nucleon numbers, respectively. 

The photonuclear cross sections a £ , (£) and o £2 (£) are Lorentzian shaped and somewhat 
sharply peaked in energy. Therefore, photon spectral functions can be taken outside of the integral of 
equation (7.46) to yield the approximate form (ref. 76) of 

CT em ~ N E\ ( £ GDr) J°£l ^ dE + N E2 ( £ GQr) £ GQrJ°£2 ^ ^2 (7.53) 
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where /^dr and £ G q R are the energies at the peaks of the E\ and E2 photonuclear cross sections, 
respectively. These integrals of photonuclear cross sections, integrated over energy, are evaluated with 
the following sum rules (ref. 76): 


r N pZp 

\°E\ dE = 6 °-T~ (7.54) 

A P 

and 

p = 0 . 22 / Z ^ (7.55) 

In equations (7.54) and (7.55), N p is the number of neutrons, Z p is the number of protons, and A p 
is the mass number of the projectile nucleus. The fractional exhaustion of the energy-weighted sum rule 
in equation (7.55) (ref. 68) is 


0.9 

(Ap> 100) 


0.6 

(40 < Ap < 100) 

(7.56) 

0.3 

(40 < Ap) 



In equation (7.53), £ GDR and £ GQR are the energies at the peaks of the E\ and E2 photonuclear 
cross sections, respectively. For the dipole term (ref. 68), 


P 

^GDR 


he m c2 ^( 
2n L 87 v 


1 + 11- 


1 + £ + 3 u 

£ 

1 +£ -I - U 



(7.57) 


with 


37 . 1/3 
M = — Ap Wi 

Q p 


(7.58) 


and 


/? - r 4 1/3 

A 0 r 0 A P 


(7.59) 


where £ - 0.0768 , Q - 17 MeV , J - 36.8 MeV , = 1.18 fm , and m* is 0.70 of the nucleon mass. 

For the quadrupole term, 


63 

: gqr - 7775 

A p 


(7.60) 


Finally, the single-proton or single-neutron removal cross sections are obtained from o 
(eq. (7.53)) by using proton- and neutron-branching ratios g. ( i = p, n). Then C 


°(0 = 8^ tm O' = porn) 

The proton-branching ratio has been parameterized by Westfall et al. (ref. 68) as 


g p = min 


— , 1.95 exp(-0.075Z ) 
L A P 


(7.61) 


(7.62) 
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where Z p is the number of protons, and the minimum value of the two quantities in brackets is to be 
taken. This parameterization is satisfactory for heavier nuclei ( Z p > 14); however, for light nuclei, the 
following branching ratios are used instead: 


For neutrons, the branching ratio is 


0.5 

(Z p <6) 


0.6 

(6<Z p < 8) 

(7.63) 
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s n = 

l ~b 

(7.64) 


Figure 3 shows the absorption cross section as a function of energy for six different projectile 
charge numbers Z p , as calculated by HZETRN for aluminum, water, and liquid hydrogen targets. The 
importance of the target size on the systematic variation of the absorption cross section with projectile 
mass is clearly seen. Small target nuclei preferentially fragmentize heavy ions, which is important to 
shield performance. (See ref. 78.) 

Among the best known fragmentation cross sections are those for carbon ion beams on carbon tar- 
gets. Measurements have been made at four energies (250A MeV in ref. 79, 600A MeV in ref. 80, and 
1.05A and 2.1 A GeV in ref. 81) and are compared in figure 4 with the results from NUCFRG2 cross- 
section calculation. The effects of the Coulomb trajectory are clearly apparent on the lighter mass frag- 
ments of lithium and beryllium which exhibit a change in slope below 100A MeV. Such Coulomb 
effects will be even more important for projectile and fragments of greater charge. Figure 5 shows the 
results of very- low-energy (1 1.7A MeV) oxygen projectiles onto a molybdenum target (ref. 82) where 
Coulomb effects are very important. The resulting charge removal cross sections seem well represented 
by NUCFRG2 calculation even at low energies. The process whereby a nucleon is exchanged between 
the oxygen projectile and the target nucleus is referred to as “exchange poles.” As shown in figure 5, 
the addition of exchange poles to the model will bring charge removal cross sections into agreement as 
can be judged by the proton exchange pole (neutron removal) contribution for proton exchange. 

Figure 6 shows three projectile-target combinations for which two groups of experimenters have 
taken measurements at nearly the same energy (1.55A GeV in ref. 83 and 1.88A GeV in ref. 84). On the 
basis of NUCFRG2 calculations, very small cross-section differences are expected to exist at these high 
energies. (See fig. 4.) Cross sections NUCFRG2 tend to favor the data of Westfall et al. (ref. 84) and lie 
10 to 50 percent above that of Cummings et al. (See ref. 83.) Perhaps the most encouraging result of the 
comparison with the data of Cummings et al. is that trends in NUCFRG2 below aluminum measure- 
ments (Z p = 13) of Westfall et al. appear in reasonable agreement with the results of Cummings et al. 
down to neon fragments (Z f - 10). The NUCFRG2 fragmentation cross sections below neon remains 
untested. 

To better quantify the comparison of results shown in figure 6, x 2 analysis is used. The test data of 
the fragmentation model in NTJCFRG2 are compared with the experimental data of Cummings et al. 
(ref. 83) and Westfall et al. (ref. 84) and shown in table I for iron projectiles on three targets. Shown in 
table I are the total % 2 values and the average x 2 contributions per degree of freedom n. Clearly, the 
data for producing aluminum fragments in the experiments of Westfall et al. show large systematic 
errors, which are nearly the sole contribution to the y} values. The elimination of the aluminum data 
from the experiments of Westfall et al. is shown in the table as the values in parentheses. By excluding 
the aluminum data, the model then shows good agreement with the data of Westfall et al. for carbon and 
copper targets. The greater discrepancy for the lead targets surely results from the simplified nuclear 
matter distribution; the use of a diffusive surface model rather than the uniform sphere model is recom- 
mended. When comparing the model with the data of Cummings et al., similar trends are shown with 


28 



target mass, but the overall agreement with the data of Cummings et al. is inferior to the agreement with 
the data of Westfall et al. as was previously noted in the discussion of figure 6. 


Table I. X 2 Analysis of Iron Fragmentation Model 


Shield + target 

NUCFRG2-Westfall a 

NUCFRG2-Cummings 

X 2 

X 2 /n 

X 2 

x 2 / n 

Fe + C 

50.2 (16.0) 

5 (1.8) 

48.3 

3.7 

Fe + Cu 

200.6 (22.9) 

20 (2.5) 

78.1 

6.0 

Fe + Pb 

177.4 (56.2) 

18 (6.2) 

83.1 

6.7 


a All values in parentheses exclude Westfall’s aluminum cross section. 


The X" analysis has been used to compare how well the results of one experimental group compare 
with the results of another. Such a comparison is shown in table II. Again, the values in parentheses 
eliminate the aluminum data contribution of Westfall et al. What is clear from table II is that the model 
represents the two sets of experimental data better than either experimental data set represents the other. 
This is the best that can be hoped for within the present experimental uncertainty. 


Table II. x 2 Analysis of Iron Fragmentation Experiments 


Shield + target 

WestfalP-Cummings 

Cummings-Westfall 3 

X 2 

X 2 /" 

X 2 

X 2 /" 

Fe + C 

85.6 (43.3) 

8.6 (4.8) 

54.6 (33.6) 

5.5 (3.7) 

Fe + Cu 

424.4 (108.4) 

42.4 (12.0) 

160.3 (69.4) 

16.0 (7.7) 

Fe + Pb 

348.8 (79.5) 

34.9 (8.8) 

143.4 (55.8) 

14.3 (6.2) 


a All values in parentheses exclude Westfall’s aluminum cross section. 


Figures 7(a), 7(b), and 7(c) show the fragmentation cross sections as calculated by NUCFRG2 as a 
function of projectile and fragment mass numbers at four different projectile energies for aluminum, 
water, and liquid hydrogen targets, respectively. It is difficult to interpret the data directly in terms of 
shield performance. However, the ion fragmentation at lower energies is clearly most effective for 
complex target nuclei as can be seen by comparing the aluminum cross sections at 25 A MeV in 
figure 7(a) with those of liquid hydrogen in figure 7(c). In contrast at high energies, the hydrogen target 
distributes the mass of the fragments more broadly than the aluminum targets as seen by comparing the 
results at 2400A MeV in figures 7(a) and 7(c). The water target displays both of these characteristics 
simultaneously. 

8. Environmental Model 

Calculations and measurements show that HZE ions in GCR particle fluxes vary inversely with 
solar activity. This dependency indicates that, at maximum solar activity, the GCR particle flux is at a 
minimum, and at minimum solar activity, the GCR particle flux reaches a maximum. This process of 
flux variation results from the shielding of the inner part of the solar system by the magnetic fields that 
are being carried away from the solar surface by the solar wind and varies periodically during an 1 1 - to 
13-year cycle. This shielding process is particularly effective in low-energy GCR particles. Calculations 
have shown that, for energies of 1004 MeV or less, the GCR particle flux variation between solar mini- 
mum and solar maximum is a factor of 10 or more. At GCR particle energies above a few hundred 
A MeV, the effect of solar modulation on GCR particle flux gradually decreases and above approxi- 
mately 10A GeV becomes negligible. 
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A numerical description of the previous solar modulation must be supplied as an input to the envi- 
ronmental model of any space radiation transport computer code including HZETRN. The environmen- 
tal model of Badhwar and O’Neill is used in HZETRN. (See ref. 42.) This model is based on fitting 
measured differential energy spectra of 1954—1989 to stationary Fokker-Planck equation solutions to 
estimate the diffusion coefficient or the equivalent deceleration parameter O in units of MV. 
Reference 42 shows that, independent of energy, this approach fits the available data within an error of 
±10 percent rms. By using the calculated diffusion coefficient, the GCR spectra of hydrogen, helium, 
carbon, oxygen, silicon, and iron are estimated. The spectra of the remaining elements were scaled to 
the previous elements by following reference 43. The environmental model data (i.e., solar maxima- 
minima) available in HZETRN are for the years 1958-1959, 1970-1971, 1981-1982, and 1989 solar 
maxima and 1965, 1977, and 1986-1987 solar minima. 

Based on the model of reference 42 used in HZETRN, figure 8 shows the predicted fluxes of 1977 
solar minimum and 1981 solar maximum as a function of energy for five different charge groups. 
Figure 9 shows the flux ratio of the 1981 solar maximum compared with the 1977 solar minimum for 
five different charge groups in the energy range of 0.1 to 10 6 MeV. As discussed previously, the two 
spectra are essentially identical at energies above 50A GeV, while at lower energies, the ratio can be 
1:10 or less depending on the ion type. 

Solar flare occurrence is correlated with solar activity and the most important events seem to occur 
during the ascending or descending phase of the solar cycle. The four solar particle events (SPE) of 
February 1956, November 1960, August 1972, and October 1989 are widely used to estimate flare- 
shielding requirements. Fitted fluence energy spectra of these four events (refs. 85-88), King’s-fitted 
spectmm of 1972, and Webber’s fit (ref. 89) are as follows: 

1. February 1956 

$p{>E) = 1.5 x 10 9 exp(- ^ 5 ^) + 3x 1° 8 ex p(~ 

2. November 1960 

<t> p (>E) = 7.6 x 10 9 exp(- + 3.9 x 10 8 exp^- E 8Q 100 ) 

3. August 1972 

4>p(>£) = 6.6 x 10 8 exp(- - 3 ^°° ) 


4. October 1989 


4 ) p (>E) = 8.65 x 1 0 10 exp 


P(E) 1 
59.261. 


5. August 1972 (King spectmm) 


V >£ ) 


7.9xl0’exp(-^) 


6. Webber spectmm with 1 00 MV rigidity 


<M> £ ) 


1.0 x 10 9 exp 


'239.1 -P(E) 
100 


where E is the energy in MeV, P{E) = Je (E + 1876) is the rigidity, and <b p ( >E) is the proton 
integral fluence in protons/cm 2 . Figure 10 shows the fitted spectra. The flare of October 1989 produced 
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the greatest number of protons less than 4 MeV, and the flare of August 1972 produced the greatest 
number of protons between 4 and 150 MeV. The February 1956 flare produced approximately one tenth 
as many protons greater than 10 MeV as the 1972 and 1989 flares, but it delivered far more protons of 
200 MeV or greater than the three other flares. The November 1960 flare exhibited intermediate 
characteristics. 

9, HZETRN Benchmarking 

In this section, the subject of HZETRN validation is addressed. Ideally, validation should be accom- 
plished with detailed transport data obtained from carefully planned and controlled experiments; unfor- 
tunately, such data are scarce. Although useful for comparative purposes, the atmospheric propagation 
measurements (ref. 36) used previously are clearly not definitive because they consist of integral flu- 
ences of as many as 10 different nuclear species combined into a single datum. Although limited quan- 
tities of HZE dosimetry measurements from manned space missions (e.g., Skylab) are available 
(ref. 90), numerous assumptions about the relationship between the dosimeter location and spacecraft 
shield thickness and geometry must be made to estimate astronaut dose with GCR codes. Because many 
of these assumptions may involve inherently great uncertainties (i.e., factors of 2 or greater), differences 
in results are difficult to attribute to the particular assumptions or approximations that may have been 
used in the analysis. Without definitive GCR transport measurements with which to compare code pre- 
dictions, other methods of validation must be considered. One such method is to compare HZETRN 
with limited available proton transport Monte Carlo results. (See ref. 87.) 

In reference 87 comparisons of the results obtained for a hypothetical problem with four different 
proton transport codes are presented. The hypothetical problem was to determine the dose, as a function 
of depth in water, resulting from a typical solar flare proton spectrum normally incident on a semi- 
infinite slab shield which is followed by the slab of water. The water was assumed to be 30 g/cm 2 thick. 
The solar proton spectrum was taken to be exponential in rigidity (Webber spectrum of section 8) with a 
characteristic rigidity of 100 MV and was normalized to 10 9 protons/cm 2 with energy >30 MeV. The 
Nucleon Transport Code (NTC) (ref. 91) is used here and compared with HZETRN. 

Figures 1 1(a) and 1 1(b) show the total dose in gray (Gy) and total dose equivalent in sievert (Sv) for 
20 g/cm shields of aluminum and iron, respectively, followed by 30 g/cm 2 of water. For both shields, 
the code correlation is good. At water depths >20 g/cm 2 , HZETRN predictions are slightly greater than 
NTC because of the 400-MeV cutoff in the NTC calculation. (See ref. 37.) 

Figures 12(a) and 12(b) show the secondary neutron dose in gray (Gy) and secondary neutron 
dose equivalent in sievert (Sv) for 20 g/cm 2 shields of aluminum and iron, respectively, followed by 
30 g/cm 2 of water. For both shields, the code correlation is only fair, which results partly from the grid 
coarseness of the Monte Carlo secondary neutron prediction runs and current limitations in the neutron 
transport database and elastic scattering propagation. (See ref. 33.) 

Other comparisons of HZETRN with analytic functions are reported in references 35, 39, and 75. 

10. HZETRN Computational Results 

In this section results from applying HZETRN to different targets and environmental periods (i.e., 
solar minimum or maximum) are presented to provide a sufficiently broad database for code testing. 
The emphasis is on information that can be extracted directly from a typical execution of HZETRN. 
Therefore, environmental periods of the 1977 solar minimum and 1981 solar maximum and target mate- 
rials of aluminum, water, and liquid hydrogen are used. Although HZETRN-calculated dose equivalent 
results are based on the International Commission on Radiological Protection ICRP-26 and ICRP-60 
(refs. 92 and 93) quality factor recommendations, only the dose equivalent results for ICRP-60 are 
shown here. The calculations include the following results for skin dose and blood-forming organ 
(BFO) dose: 
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1 . Flux contributions from carbon, oxygen, calcium, and iron ions at different energies 

2. Total and individual ions, dose and dose equivalent contributions for skin (0-cm surface dose), 
and BFO (5-cm-depth dose) at different depths 

3. Flux contributions, total dose and dose equivalent contributions for skin (0 cm), and BFO (5 cm) 
at different linear energy transfers (LET) 

Major shortcomings of the HZETRN calculations are as follows: 

1. All secondary particles from HZE interactions are presently assumed to be produced with a 
velocity equal to that of the incident particle; this is conservative for neutrons produced in HZE 
particle fragmentations 

2. Meson contributions to the propagating radiation fields are neglected 

3. Target fragments in HZE reactions are neglected 

4. Angular dependence, especially in neutron propagation is neglected 

These data are not all conservative and probably account for the 15- to 30-percent underestimation 
of the exposure. As discussed in references 94 and 95, the main sources of uncertainty are the input 
nuclear fragmentation model and the incident GCR spectrum. Taken together, they could easily impose 
an uncertainty factor of 2 or more in the astronaut risk estimate. This is especially true for track struc- 
ture dependent biological models. (See refs. 78 and 96.) 

Figures 13-16 show the flux of carbon, oxygen, calcium, and iron ions in aluminum targets for the 
1977 and 1981 solar periods at depths of 0, 10, 20, and 30 g/cm 2 in the energy range of 10 2 A to 
10 5 A MeV. Similarly, figures 17-24 show the flux for water and liquid hydrogen targets. In all figures 
the effect of solar minimum and solar maximum periods are demonstrated through a reduction in the 
magnitude of the particular ion flux during solar maximum. Figures 13-24 demonstrate that lighter 
target materials such as liquid hydrogen have better attenuation characteristics than heavier target 
materials such as aluminum for all target depths. 

Figure 25 shows skin dose, Gy and dose equivalent, Sv, when using ICRP-60 quality factors 
in aluminum, water, and liquid hydrogen targets for the 1977 and 1981, 1982 solar periods at depths of 
0 to 30 g/cm 2 . Only a modest gain in protection is displayed by aluminum shielding. A water shield 
shows substantially better attenuation characteristics. Liquid hydrogen has excellent shielding charac- 
teristics. These results are demonstrated more clearly in figure 26, which shows skin dose, Gy, and dose 
equivalent, Sv, when using ICRP-60 quality factors in aluminum, water, and liquid hydrogen targets for 
the 1977 and 1981 solar periods at depths of 0 to 30 g/cm 2 . In figure 26, the effect of the atomic and 
nuclear properties of the target is demonstrated by increases in skin doses for larger target nuclei. 

Figures 27-29 show skin dose, Gy, and figures 30-32 show dose equivalent, Sv, when using 
ICRP-60 quality factors for six different charge groups in aluminum, water, and liquid hydrogen targets 
for the 1977 and 1981 solar periods at depths of 0 to 30 g/cm 2 . The better attenuation properties of the 
liquid hydrogen target versus aluminum target at all target depths are shown for the given charge groups 
for both solar epochs. 

Figures 33-40 show the attenuation behavior of GCR in aluminum, water, and liquid hydrogen at 
the 5-cm-BFO depth. The dosimetric characteristics of different target materials as a function of target 
depth are similar to those in figures 13-32, except here an additional attenuation exists in all data 
because of an additional 5 g/cm 2 of water. 

Figures 41-49 show skin flux and dosimetric characteristics of GCR components in aluminum, 
water, and liquid hydrogen targets for the 1977 and 1981 solar periods for the LET range of 1 to 
10 5 MeV/cm of water at depths of 0, 10, 20, and 30 g/cm 2 . The presence of the primary components of 
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GCR from proton to iron is clearly demonstrated as inflections in all figures. The attenuation character- 
istics of lighter targets such as liquid hydrogen versus heavier targets such as aluminum are shown at all 
depths and LET ranges. 

Figures 50-58 show 5-cm-BFO flux and dosimetric characteristics of GCR components in alumi- 
num, water, and liquid hydrogen targets for the 1977 and 1981 solar periods for the LET range of 1 to 
10 MeV/cm of water at depths of 0, 10, 20, and 30 g/cm 2 . As before, the dosimetric characteristics of 
different target materials as a function of target depth are similar to those in figures 41-49, except here 
an additional attenuation exists in all figures because of an additional 5 g/cm 2 of water. As in 
figures 41-49, the presence of the primary components of GCR from proton to iron is clearly demon- 
strated as inflections, and the attenuation characteristics of lighter targets such as liquid hydrogen versus 
heavier targets such as aluminum are shown at all depths and LET ranges. 

The use of the HZETRN code and the control of the options to generate a typical output are 
discussed in appendix A. A complete sample output with the options of appendix A is provided in 
appendix B. 

11. Concluding Remarks 

The study of free-space radiations is essential for estimation of the risk factors for long-duration 
manned space flights. The estimates of galactic cosmic ray (GCR) ion flux inside a space module can be 
made by calculating quantities such as energy and linear energy transfer (LET) spectra. The estimation 
of doses and LET spectra by numerical predictive techniques, such as HZETRN, is essential for any 
long-duration manned space mission. 

This paper reports the progress in computer code development for space radiation studies and shield 
design for future NASA space programs; the HZETRN code is a state-of-the-art fast computational tool 
available for a design engineer to obtain answers to some of the radiation questions that arise in plan- 
ning any mission. However, major uncertainties in nuclear cross sections, environmental models, and 
astronaut risk affect the overall accuracy of the predictions of any analytical-computational technique. 
These uncertainties have a major impact on the proposed shield design for any mission and the subse- 
quent mission cost. Much work remains to accurately resolve the problems with nuclear cross-section 
calculations, environmental model development, and risk estimate methods. 


NASA Langley Research Center 
Hampton, VA 23681-0001 
February 14, 1995 
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Appendix A 

HZETRN Description 

Al. Program Structure 

The main program and each subroutine or function module begins with a brief description of its pur- 
pose. The program has 8358 lines, is written in FORTRAN 77, and was developed on a VAX 4000 
scaler minicomputer under the VMS operating system; in its current version, the program requires 
463 blocks (237056 bytes) of storage space. However, the program is essentially platform and operating 
system independent. The complete computer code consists of a main program HZETRN, 63 subroutines, 
and 37 function modules. Table AI provides a list of the subroutines and function modules in the order 
in which they appear in the program. The overall structure of HZETRN is shown in figure Al, which 
includes input requirements and output options. 


Table AI. HZETRN Program Modules 


No. 

Program 

module 

Type 3 

No. 

Program 

module 

Type 3 

No. 

Program 

module 

Type 3 

1 

MATTER 

s 

35 

DENSEFF 

F I 

69 

FB 

F 

2 

DMETRIC 

s 

36 

ATOPA 

s 

70 

SLOPE 

F 

3 

PRPGT 

s 

37 

ACOEF 

F 

71 

ELSEC 

F 

4 

OD 

F 

38 

CCOEF 

F 

72 

XTOT 

F 

5 

TS 

F 

39 

ATOPN 

s 

73 

XTOTLM 

F 

6 

F2 

S 

40 

BUILD 

F 

74 

XSEC 

F 

7 

FRAG 

F 

41 

SNF 

F 

75 

ABSEC 

F 

8 

XSECFRG 

F 

42 

TEXP 

F 

76 

LZEVAP 

S 

9 

YIELDX 

S 

43 

TSQR 

F 

77 

ASIGM 

S 

10 

YIELDN 

S 

44 

RMAT 

F 

78 

XTEST 

S 

11 

YIELDA 

S 

45 

SIONA 

F 

79 

SPLCN 

S 

12 

YIELDT 

S 

46 

ZEFF 

F 

80 

SPLIN 

S 

13 

PXN 

S 

47 

RTIS 

F 

81 

PRPLI 

S 

14 

YIELDl 

S 

48 

MGAUSS 

S 

82 

FD 

s 

15 

YIELD2 

S 

49 

I UNI 

s* 

83 

LITEST 

s 

16 

YIELD3 

s 

50 

IBI 

s* 

84 

TBARLI 

s 

17 

YIELD4 

s 

51 

QXZ114 

F* 

85 

LIFRAG 

s 

18 

NAND 

F 

52 

QXZ115 

s* 

86 

RADIUS 

F 

19 

YIELDH 

S 

53 

QXZ116 

s* 

87 

XLISEC 

F 

20 

YIELDLI 

S 

54 

QXZ117 

s* 

88 

LKOSPC 

S 

21 

YIELDEM 

F 

55 

SUTPAR 

s* 

89 

E3MAX 

F 

22 

BESSEL 

F 

56 

QXZ060 

s* 

90 

WDKO 

S 

23 

GEODA 

S 

57 

QXZ061 

s* 

91 

LQESPC 

s 

24 

BSEACH 

s 

58 

QXZ062 

s* 

92 

WQE 

F 

25 

DELTA 

s 

59 

QXZ063 

F* 

93 

DERF 

F 

26 

FPOFB 

s 

60 

QXZlll 

s* 

94 

GCRFLUX 

S 

27 

LIMIT 

s 

61 

QXZ112 

s* 

95 

FBERT 

s 

28 

GEOFR 

s 

62 

QXZ119 

s* 

96 

BERT INI 

s 

29 

RADLI 

F 

63 

PHTOINT 

s 

97 

RANFT 

s 

30 

EOFS 

S 

64 

PHI 

F 

98 

FDNP 

s 

31 

CROSSO 

s 

65 

PHTOD 

S 

99 

Q60 

F 

32 

QOFS 

F 

66 

ELSPEC 

S 

100 

SFBERT 

s 

33 

ATOE 

S 

67 

XPP 

F 




34 

ATOPP 

s 

68 

XNP 

F 





a S — subroutine module. 

F — function module. 

S* — subroutine module taken from LaRC-NASA mathematical library (N2-3D). 
F* — function module taken from LaRC-NASA mathematical library (N2-3D). 
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A2. HZETRN Input Requirements 

The complete program is composed of the HZETRN source code (HZEALH20 . FOR in the COSMIC 
library) and three input data files. In the current version, HZEALH20 . FOR requires the following input 
files: 

1. ATOMIC . DAT 

2. NUCLEAR . DAT 

3. JSCCOSMC . DAT 

The following is a brief description of each file and how to set up the input flags if any are required. 

1. FORTRAN code HZEALH20 . FOR is set up to transport galactic cosmic ray (GCR) particles in 
free space (geomagnetic cutoffs are ignored) through a given thickness of the aluminum shield followed 
by a given depth of water. Besides the setting of array dimensions for the energy grid points and isotope 
fragment number, which are defined in the PARAMETER statement throughout the code as 11=45 and 
IJ=59, respectively, the only changes to the code required for any run are in the MAIN module and 
subroutine MATTER. In addition the user has the option of changing the isotope table that is used to rep- 
resent the transported particle fields of the GCR particles in the shield or target material in function TS. 

a. In the MAIN module, the following two integer variables must be defined: 

(1) IPRGCR, which specifies the extent (1, 2, or 3) of printed information, with 2 being a good 
choice for most calculations 

(2) IYEAR, which specifies the year in the solar cycle to be used (solar maxima for the years 
1958, 1970, 1981, and 1989 and solar minima for the years 1965, 1977, and 1986) 

b. In subroutine MATTER, the following arrays and variables must be defined: 

(1) Dimension of array XAL, which specifies the depth in the aluminum shield where dosimetric 
quantities are to be calculated and printed. The depth within the aluminum shield is deter- 
mined by summing the elements of array XAL. For example, if XAL has a dimension of 3 and 
is defined in the data statement as DATA XAL/ 5 , 5 , 10/, the transport calculation results 
will be printed for the following depths of the shield material: 

(a) Depth 1 : XAL ( 1 ) = 5 g/cm 2 

(b) Depth 2: XAL ( 1 ) +XAL (2) =5 + 5 = 10 g/cm 2 

(c) Depth 3: XAL ( 1 ) +XAL ( 2 ) +XAL ( 3 ) =5 + 5 + 10 = 20 g/cm 2 

(2) Dimension of array XWAT. The discussion in paragraph l.b.(l) also applies to the array XWAT 
for water. 

Note that for the shield and target material the particle fields are transported and calculated by 
default in steps of H=1 g/cm 2 throughout the program, and the discussion in para- 
graphs 1 .b.(l) and 1 .b.(2) refers only to the depths for which data are to be printed. 

(3) Variables NTOT1 and NTOT2 must be equal to the dimensions of XAL and XWAT arrays, 
respectively. For example, if XAL has dimension 3 and XWAT has dimension 5, then NTOT1 
and NTOT2 are defined by 


DATA NTOT1 , NTOT2 /3 , 5/ 
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(4) The following physical properties of the shield and target material are included: NAT is the 
number of elements in the shield or target — NAT = 1 for aluminum, NAT =2 for water, DN is the 
overall density of the shield or target in g/cm 3 , ATRG is the atomic mass number of each ele- 
ment of the shield or target, ZTRG is the atomic charge number of each element of the shield 
or target, and DENSTRG is the number of atoms/g of each element of the shield or target. 
DENSTRG is calculated from 

DENSTRG (I)=[N0*R(I) ] /SUM [R ( I) *ATRG ( I ) ] 

where NO is the Avogadro number taken as 6 . 023E23 atoms/gram-atom, R ( I ) is the num- 
ber of atoms/molecule, and ATRG ( I ) is the mass in gram/gram-atom of the Ith element in 
the shield or target. For a simple material like aluminum, DENSTRG becomes 

DENSTRG ( 1) = ( 6 . 023 E2 3*1) / (1*27.0) =2. 23 E2 2 (for Al) 

For a compound material like water, DENSTRG becomes 

DENSTRG ( 1 ) = (6 . 023E23*2 ) / (2*1 . 008+1*16 . 0) = 6 . 68E22 (for H) 

DENSTRG (2 ) = (6. 023 E23*l)/(2*1. 008 + 1*16. 0)=3 . 3 4E22 (for O) 

Note that only 3 significant figures are retained for compatibility with the generated data files. 
Also, note that arrays ATRG, ZTRG, and DENSTRG currently have a dimensional limit of 5. 
For a more complex material with NAT>5, the dimensions of ATRG, ZTRG, and DENSTRG 
must be adjusted appropriately. 

(5) The transport of chosen particles in the shield or target is represented by a user-provided iso- 
tope table that is defined in function TS. The arrays APRO and ZPRO in common block PROJ 
contain the mass and charge numbers of transported isotopes. The isotope table can be 
changed by choosing a different table or, if more or less than 59 isotopes are needed, by 
adjusting variable IJ=59 in the PARAMETER statement throughout the codes to reflect a 
larger or smaller isotope table. 

2. ATOMIC . DAT is the energy, range, and stopping-power database for water and aluminum. For 
any other material, the source code will automatically generate a file called NEWRMAT . DAT during exe- 
cution. This file should be appended to the ATOMIC . DAT master data file for future transport runs. 
Note that DENSTRG format is E12 . 3, which permits only three significant figures. 

3. NUCLEAR . DAT is the absorption and fragmentation cross-section database for water and alumi- 
num. For any other material, the source code will automatically generate a file called NEWNUC5 9 . DAT 
during execution. This file should be appended to the NUCLEAR . DAT master data file for future trans- 
port runs. Note that DENSTRG format is E12 . 3, which permits only three significant figures. 

4. JSCCOSMC . DAT is the environmental model database taken from references 42 and 43 for the 
solar maxima for the years 1958, 1970, 1981, and 1989 and solar minima for the years 1965, 1977, and 
1986. 

A3. Compiler Dependencies 

Besides the usual compiler-dependent OPEN and REAL* 8 statements, which must be checked by 
the HZEALH20.FOR code user, the source code has one compiler-dependent feature in subroutine 
QXZ0 61, which must be modified if this code is to be run on a non-VAX platform. In QXZ0 61, 
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variables THIRD and SIXTH are defined in OCTAL format and must be modified to conform with a 
non- VAX compiler octal number format. Besides these trivial changes, the rest of the source code 
should run on any platform with a FORTRAN compiler. 

A4. Program Execution 

To provide the users of HZEALH20 . FOR with a detailed description of a typical run, the following 
input conditions are used to generate the sample output given in appendix B: 

1. IPRGCR=2 

2. IYEAR=4 (1977 solar minimum) 

3. XAL ( 3 ) = 5 , 5 , 1 0 (aluminum depths of 5, 10, and 20 g/cm 2 ) 

4. XWAT (5) =1,1, 1,1,1 (water depths of 1, 2, 3, 4, and 5 g/cm 2 ) 

The sample output is divided into 36 blocks (A Z and AA JJ). Each block includes a brief 

description of the flow of the program and how to interpret the output. 

Block A prints the values of the physical properties (mass, charge, density, etc.) of the target mate- 
rial, which in almost all cases is water. These quantities are used in the function RTIS. In RTIS the 
input data file ATOMIC . DAT is accessed to see if it contains the energy, range, and stopping power of 
the chosen target material. This is accomplished by comparing the supplied physical properties of the 
target with those available in ATOMIC . DAT. If the physical properties are the same, then the energy, 
range, and stopping-power arrays are loaded from ATOMIC . DAT. If ATOMIC . DAT does not have the 
information for the target of interest, then a host of other subroutines and functions are called to gener- 
ate the appropriate energy, range, and stopping-power arrays. This new information is written into the 
file NEWRTIS . DAT. The user has the responsibility to append NEWRTIS . DAT to ATOMIC . DAT after 
the program execution for future use. The user is then informed that all appropriate arrays in RTIS are 
initialized, and the message RTIS IS INITIALIZED is printed. 

Next, the user is informed that the 1977 solar minimum is chosen as the environmental model. 
The program searches the input data file JSCCOSMC.DAT to load and scale the flux arrays in 
HZEALH20 . FOR for later transport calculations. 

Then, the program reads the physical properties (same quantities as previously read in RTIS) of the 
shield material from subroutine MATTER. The subroutine MATTER, as discussed previously, is the only 
module that the user has to modify to accommodate new shield and target materials. Once the shield 
material information is read from MATTER, the information concerning the transported particle field 
absorption and fragmentation cross sections for the shield material is read from the input data file 
NUCLEAR . DAT. This is accomplished by calling function TS. In TS, as in RTIS, the physical proper- 
ties of the aluminum shield material are compared with those available in NUCLEAR . DAT. If the phys- 
ical properties are the same, then the absorption and fragmentation cross-section arrays are loaded from 
NUCLEAR . DAT. If NUCLEAR . DAT does not have the information for the shield material of interest, 
then a host of other subroutines and functions are called to generate the appropriate cross-section arrays. 
This information is written into the file NEWNUC59 . DAT. As in RTIS, the user has the responsibility to 
append NEWNUC59 . DAT to NUCLEAR . DAT for future runs. The user is then informed that all appro- 
priate cross-section arrays are initialized and the message MATERIAL FOR TS 1 FOUND ON 
NUCLEAR . DAT FILE is printed. 

The physical properties of the shield material are used in function RMAT to load energy, range, and 
stopping-power arrays for the shield material from input data file ATOMIC . DAT. RMAT is identical to 
RTIS, except that it is designed to calculate or load the energy, range, and stopping-power arrays for 
shield material rather than for target material. As in RTIS, if the ATOMIC . DAT does not have the infor- 
mation for the shield material, a new file called NEWRMAT . DAT is generated which will have to be 
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appended to ATOMIC . DAT for future runs. Finally, when all appropriate arrays in RMAT are initialized, 
the message RMAT IS INITIALIZED is printed. 

Block B contains the depth values of the shield material in g/cm 2 at which transported quantities are 
printed. The first depth value corresponds to the free-space boundary or X=0 g/cm. The quantities for 
the transported particles are atomic mass, atomic charge, integral flux for a given isotope in particles/ 
cm 2 -yr, dose in cGy, and dose equivalent in cSv according to the International Commission on Radio- 
logical Protection ICRP-26 and ICRP-60 guidelines are printed for charge Z=0 to Z=2 8. Note that dose 
and dose equivalent can be easily converted to engineering units by the following: 

1 rad = 1 centiGray (cGy) = 0.01 Gy 
1 rem = 1 centiSievert (cSv) = 0.01 Sv 

Then, the minimum and maximum energy, A MeV, and range, g/cm 2 , that were used to calculate the 
previous flux quantities are printed. 

Block C prints the cumulative integral flux and dosimetric values at the depths of the shield material 
that correspond to the same depths as block B. The cumulative transported quantities are grouped into 
six charge groups of Z=0, Z=l, Z=2, 3<Z<10, 11<Z<20, and 21<Z<28. The units for integral flux, 
dose, and dose equivalents are the same as for block B. 

Then, the total cumulative integral flux and dosimetric values for all six charge groups are printed. 

Block D contains the depth values of shield material at which transport quantities are calculated and 
stored. The depth increment is set by default to H=1 g/cm 2 . The printing of depth values stops when the 
first depth of the shield material as defined in subroutine MATTER is reached. 

Blocks E, H, and K are similar to block B. 

Blocks F, I, and L are similar to block C. 

Blocks G and J are similar to block D. At this stage, calculations for the aluminum shield material 
are complete and the program begins to read information for the target material. 

Block M prints the values of the physical properties (mass, charge, density, etc.) of the target mate- 
rial, as defined in subroutine MATTER. These quantities are used in function RMAT. As described in the 
second half of block A, these physical properties are compared with those available in the input data file 
ATOMIC . DAT, and the appropriate arrays for energy, range, and stopping power are either loaded from 
ATOMIC . DAT or calculated and stored in NEWRMAT . DAT. As in block A, the user has the responsibil- 
ity to append NEWRMAT . DAT to ATOMIC . DAT after the program execution for future use. 

Blocks N, Q, S, U, W, and Y are similar to block B. 

Blocks 0, R, T, V, X, and Z are similar to block C. 

Block P is similar to block D. At this stage, calculations for the water target material are complete 
and the program begins the stopping-power calculations for the shield and target materials. 

Block AA is the start of the printed transported quantities as a function of stopping power S=LET. It 
contains the depth values of the shield material in g/cm 2 , at which the transported quantities are calcu- 
lated. The first depth value corresponds to the free-space boundary at X=0 g/cm 2 . The printed quantities 
for the transported particles are LET in MeV/(g/cm 2 ), integral flux at a given LET in particles/cm 2 yr, 
dose in cGy, and dose equivalent in cSv to ICRP-26 and ICRP-60 guidelines. 

Blocks BB, CC, and DD are similar to block AA. The printed quantities in blocks AA, . . . , DD are 
LET-related quantities in the shield material. 

Block EE prints the transported quantities as a function of LET for the target material at a target 
depth of X=0 g/cm 2 . 

Blocks FF, GG, HH, II, and JJ are similar to block EE. 
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Environment 
Energy grid 1 
Shield material 
Shield thickness 


Energy loss database 
Cross-section database 
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Quality factor specification 
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Alternate risk estimates 
LET spectra 


1 Designates semiflexible inputs. 

Figure A1 . Computational structure of HZETRN. 
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Appendix B 

HZETRN Computer Code Sample Output 
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EMIN . EMAX = O . 1 0000 E -01 0 . 50000 E +05 
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FOR Z = 0 FNUT.DNUTt HNUT (26) .HNUT(60)= 0.143E + 09 0.105E + 00 0.210E + 01 0.210E + 01 



FOR Z* 1 FPDT. DPDT. HPDT (26) . HPDT< 60) * 0.201E+09 0.142E+02 0.210E+02 0.209E + 02 
FOR Z*2 FALP* DALP. HALF (26) # HALP<60) ■ 0.110E+08 0.264E+01 0.793E+01 0.928E+01 

FOR 3< *Z< * 10 FL IT. DLI T. HL 1 T (26) . HL ! T < 60 ) * 0.836E+06 0.173E+01 0.908E+01 0.896E + 01 

FOR 11< = Z020 FMED. DMED. HMED (26) . HMED ( 60 ) * 0. 963E+05 0.814E+00 0.952E+01 0.128E+02 

FOR 21 < = Z< *26 FHVY * DHVY . HHVY (26) » HHVY ( 60 > = 0.177E+05 0.476E+00 0.912E + 01 0. 107E + 02 
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RELEVANT QUANTITIES AT X= 0.000 FOLLOW 
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Figure 1 . Transport of particles through spherical region. 
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(a) Range and stopping power in aluminum. 
Figure 2. Linear energy transfer (LET). 
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(a) Aluminum target. (b) Water target. 
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(c) Liquid hydrogen target. 
Figure 3. Absorption cross section cs 
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Figure 4. Charge-removal cross sections for carbon on carbon. 
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Figure 5. Charge-removal cross sections for 1 1 .7A -MeV oxygen projectile on molybdenum target. 
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(a) Carbon target. 



(b) Copper target. 


Figure 6. Charge-removal cross sections for 1 ,8&4-GeV iron projectile. 
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(c) Lead target. 


Figure 6. Concluded. 
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(a) Four energy levels in aluminum. 
Figure 7. Fragmentation cross sections. 
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(b) Four energy levels in water. 
Figure 7. Continued. 
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(c) Four energy levels in liquid hydrogen. 
Figure 7. Concluded. 
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Figure 9. GCR ion flux ratio of 1981 solar maximum to 1977 solar minimum. 
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Figure 10. Integral fluence spectra of important solar particle events. 
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(a) 1 977 solar minimum. 
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(b) 1981 solar maximum. 

Figure 14. GCR energy flux for oxygen ion in aluminum. 




(b) 1981 solar maximum. 


Figure 15. GCR energy flux for calcium ion in aluminum. 
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(a) 1 977 solar minimum. 
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(b) 1 98 1 solar maximum. 

Figure 18. GCR energy flux for oxygen ion in water. 
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Figure 22. GCR energy flux for oxygen ion in liquid hydrogen. 
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(a) 1977 solar minimum. 



Energy, A MeV 


(b) 1981 solar maximum. 

Figure 23. GCR energy flux for calcium ion in liquid hydrogen 
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Figure 25. GCR skin dose and dose equivalent. 
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Figure 25. Continued. 
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(c) In liquid hydrogen. 
Figure 25. Concluded. 
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(a) 1977 solar minimum. 
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(b) 1981 solar maximum. 

Figure 28. GCR skin dose for six charge Z groups in water. 
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(a) 1977 solar minimum. 
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(b) 1981 solar maximum. 

Figure 30. GCR skin dose equivalent for six charge Z groups in aluminum. 
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Figure 32. GCR skin dose equivalent for six charge Z groups in liquid hydrogen. 







Figure 33. Continued. 
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(c) In liquid hydrogen. 


Figure 33. Concluded. 





Figure 34. GCR 5-cm-BFO dose and dose equivalent in aluminum, water, and liquid hydrogen. 
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(a) 1977 solar minimum. 
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(b) 1981 solar maximum. 

Figure 37. GCR 5-cm-BFO dose for six charge Z groups in liquid hydrogen 
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(b) 1981 solar maximum. 

Figure 39. GCR 5-cm-BFO dose equivalent for six charge Z groups in water 
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(a) 1977 solar minimum. 
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(b) 1981 solar maximum. 

Figure 40. GCR 5-cm-BFO dose equivalent for six charge Z groups in liquid hydrogen. 
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(a) 1 977 solar minimum. 



(b) 1981 solar maximum. 


Figure 42. GCR skin integral flux LET spectra in water. 
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(a) 1 977 solar minimum. 
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(b) 1981 solar maximum. 

Figure 50. GCR 5-cm-BFO integral flux LET spectra in aluminum. 
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(b) 1981 solar maximum. 



(a) 1977 solar minimum. 
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(b) 1981 solar maximum. 


Figure 52. GCR 5 cm-BFO integral flux LET spectra in liquid hydrogen 
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Figure 54. GCR 5-cm-BFO dose LET spectra in water. 
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